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A FORTRAN PROGRAM FOR CALCULATING 
THREE-DIMENSIONAL, INVISCID, ROTATIONAL 
FLOWS WITH SHOCK WAVES IN AXIAL 
COMPRESSOR BLADE ROWS - USER’S MANUAL 


by 


William T. Thompkins, Jr. 


Gas Turbine and Plasma Dynamics Laboratory 
Department o£ Aeronautics and Astronautics 
Massachusetts Institute of Technology 
Cambridge, Massachusetts 


SUMMARY 

A FORTRAN-IV computer program has been developed for the calcula- 
tion of the inviscid transonic/supersonic flow field in a fully three-dimen- 
sional blade passage of an axial compressor rotor or stator. Rotors may have 
dampers (part-span shrouds). MacCormack’s explicTit time— marching method is 
used to solve the unsteady Euler equations on a finite difference mesh. This 
technique captures shocks and smears them over several grid points. Input 
quantities are blade row geometry, operating conditions and thermodynamic 
quantities. Output quantities are three velocity components, density and 
internal energy at each mesh point. Other flow quantities are calculated 
from these variables. A short graphics package is included with the code, 
and may be used to display the finite difference grid, blade geometry and 
static pressure contour plots on blade-to-blade calculation surfaces or blade 
suction and pressure surfaces. 

Flows in four transonic compressor rotors have been analyzed and 
compared with exit flow field measurements and intra-blade static density 
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measurements obtained with a gas fluorescence technique. These comparisons 
have generally shown that the computed flow fields accurately model the 
experimentally determined passage shock positions and overall aerodynamic 
performance. 

The computer code was developed and generally run on a large mini- 
computer system, a Digital Equipment Corporation PDP-11/70, with run times of 
two to three days. The code has also been run on several main-frame com- 
puters (IBM 3033, IBM 360/678, UNIVAC 1110, CDC 7600 and a CRAY-1). Typical 
run times on an IBM 3033 have been found to be 5-10 hours. 
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INTRODUCTION 

For many analysis problems in turbomachinery, an assumption of 
inviscid flow provides sufficiently accurate results for design or develop- 
ment tasks. This situation often arises in calculation of design point per- 
formances for high speed axial compressor blade rows, even though these flows 
may exhibit complex interactions of inviscid and viscous phenomena. An 
accurate inviscid calculation is of great benefit for compressor design since 
viscous and shock losses can both be reduced if the inviscid flow can be pre- 
dicted. The design point performance of high speed axial compressors could 
be greatly improved through the use of an accurate inviscid flow solution, 
without requiring a fully viscous flow solution. 

Classical analytical solutions for these flows have not developed 
either because of the flow character — mixed subsonic-supersonic flow with 
strong shock waves — or because of the intrinsic nonlinear, three-dimensional 
features of transonic flows. For such flows pure numerical procedures can 
usually provide quick, accurate solutions with reasonable cost. The numeri- 
cal procedures and techniques to be described have been specialized for solu- 
tions of inviscid flow in high pressure ratio, high tip-speed axial 
compressor rotors , even though the techniques employed are of much wider 
generality. 

The compressor rotors to be studied are assumed to be isolated 
blade rows which are completely enclosed by hub and tip casings . Since the 
fluid is assumed to be inviscid and each blade in a row is assumed to be 
identical, the flow field about each blade may reasonably be considered to be 
identical. This assumption allows the physical domain of interest to be 
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reduced to that bounded by a pair of blades and the extension of their mean 
camber lines upstream and downstream, as Illustrated In Figure 1. This 
region Is assumed to extend to upstream and downstream Infinity, and the flow 
Is assumed to be periodic blade-to-blade. The blade rows of interest have 
several common features: their internal flows have mixed subsonlc- 
transonic-supersonlc relative Mach numbers, a range of 0.5 to as high as 2.0, 
and they attempt to use moderately strong shock waves as an efficient method 
to transfer energy from the rotating machinery to the fluid flow; their 
Internal flow passages are complex three-dimensional shapes in which natural 
bounding surfaces rarely join orthogonally. These characteristics rather 
severely restrict the present choice of numerical solution schemes to time- 
accurate, finite difference solutions of the three-dimensional Euler 
Equations (continuity, inviscid momentum and inviscid energy equation with no 
heat conduction). 

The numerical scheme selected to Integrate the equations of motion 
1 2 

is MacCormack’s method. * This scheme is an explicit, time-accurate, con- 
ditionally stable method of second order accuracy with good shock capturing 
properties. Shock waves are resolved as regions of high gradients spread 
over about 5 mesh points in the streamwlse direction. In a complex flow 
the existence and location of shock waves need not be anticipated but 
develop naturally as the solution proceeds . The penalty paid for this con- 
venience is a loss in spatial resolution of shock waves and some inaccuracy 

in shock jump conditions. Although shock-fitting schemes are an area of 
3 4 

current research, * methods suitable for three-dimensional flows are not yet 


available . 
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In most situations, the time accurate nature of MacCormack’s method 
is of little importance since only the steady state solution is of interest. 
Here the integration method is only a convenient iteration scheme to move 
from a rather arbitrary initial condition to the final steady state solution. 
Commonly used initial conditions are an old steady state solution or a quasi- 
meridional start-up procedure provided vd.th the code. 

Finite difference methods are useful for complex equation sets such 
as the Euler or Navier-Stokes equations, but can be effectively used only in 
simple geometric regions. Using finite difference methods for complex 
geometries requires adopting coordinate systems or coordinate mappings which 
transform the physical domain into a computational domain of simple shape. 

In this report, a set of simple analytical stretching functions is used to 
map the bounding surfaces Into a square computational domain. These trans- 
forms accommodate any hub or tip casing shape and almost any blade shape. 

Best results are obtained for thin leading and trailing edges as are most 
often encountered in compressor blade rows. Further work is continuing on 
more general coordinate mappings to provide more accurate solutions for any 
blade geometry. 

The computer codes described in this report represent a radical 
departure in philosophy of large-scale computational projects in the choice 
of computer systems. These codes were developed and production runs made on 
a dedicated minicomputer, a Digital Equipment Corporation PDP-11/70, rather 
than a large main frame computer. These codes have also been run on an IBM 
360/67, an IBM 3033, a UNIVAC 1110, a CDC 7600, and a CRAY-1. It was found 
during development that, at least in the author's opinion, the dedicated 



6 


minicomputer was a superior code development machine to any main frame com- 
puter. The cost of production code running was also lower than on any 
available main frame computer, with a total cost of about 300 dollars per 
solution. 
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METHOD OF ANALYSIS 


Flow Equations 

For flow calculations in turbomachinery blade passages, a cylindri- 
cal coordinate system is a natural choice since considerable cylindrical sym- 
metry exists, and was selected as a base coordinate system. This system is a 
right-handed one (r, 9, z) which has the positive z coordinate pointing in 
the axial flow direction. Blade rows will be assumed to rotate in the posi- 
tive 9 direction, clockwise when looking downstream. 

A convenient set of three-dimensional inviscid flow equations 
expressed in cylindrical coordinates is found in MacCormack.^ These equations 
are in weak conservation law form and may be easily expressed in matrix form 
as: 
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= p(e+u /2), and u“=u^ +Ug +u^ . In addition the working fluid will be 

assumed to be a perfect gas with constant thermodynamic properties, Y> R» 

C and C , The equation of state may be written as 
p V 

p = pe(Y - 1) (2) 

In order to apply the boundary condition of no flow through the blade 
surface at a constant spatial location a common independent variable trans- 
formation is introduced* 

e’ = 0 - wt (3) 

The new independent variable is 0', and u) is the rotational speed of 
the blade row* Using this transform the flow equations become: 


at’’ 3 0’’ ar' 30' 32' 


(4) 
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where G = G* - ujU 


pUg “ (urp 

pu u« - o)rpu 
r 0 r 

2 

pu^ + P “ wrpu^ 
pu^ue - wrpu^ 
Uq(E^+P) - (orEj. 


and the prime superscripts have been dropped for convenience. 

These equations are non-dimensionalized using a reference length, 

L, a reference velocity, a^, and a reference density, p^. The reference 
quantities are in principle arbitrary, but the reference velocity and density 
have been selected to be the stagnation speed of sound and density on the 
inlet tip casing streamline. The reference length remains arbitrary but is 
conventionally selected as the inlet tip casing radius. The new non~ 
dimensional variables (primes indicate dimensionless quantities) arei 


length 

r' = r/L , z' = z/L 

(6) 

velocity 

u ' = u/a^ 

(7) 

density 

P ' = PQ 

(8) 

pressure 

p’ = p/Pq^^ 

(9) 

energy 

e* = e/a^ 
0 

(10) 

rotational speed = wL/a^ 

(11) 


The reference quantities T^ and are not independent and are 


determined from: 
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= a^/ yR 
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(12) 

= 1/y 

(13) 


\^en these quantities are selected, the flow equations in terms of non- 
dimensional variables are identical in form to those in terras of dimensional 
variables, equation (5). The prime superscripts may then be dropped without 
confusion and only nondimensional quantities will be referred to in the 
remaining sections of this report. The Euler flow equations become 
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Transformed Equations 

As is typical in three-dimensional flow solutions, the geometric 
domain in physical space is sufficiently complex to render a direct finite 
difference solution in physical space impractical. The approach adopted here 
to solve the finite difference grid allocation problem is to introduce a set 
of independent variable transformations or coordinate stretchings which map a 
single compressor blade passage into a rectangular parallelepiped. A general- 
ized mapping is indicated by 


with 


(r, e, z) > (R, 0, X) 


(15) 


R * R(r, 9, z), 0 - 0(r, 0, z), and X = X(r, 0, z) . 

Introducing a coordinate transform modifies the original partial 
differential equations. Using the chain rule for derivatives Equation (14) 
becomes 


3U 3F 9R 3F 30 3F 3X 9G 3R 3G 30 3G 3X 

3t 3R 3r 30 3r 3X 3z 3R 30 30 39 3X 30 


( 16 ) 


3R 3z 30 3z 3X 3z “ ^ 


This equation may be rearranged as: 
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3F 90 3G 90 9H 30 

90 9r 30 96 90 9z 


(17) 


3F 3X 3G 9X 9H 9X 

3X 9r 3X 30 9X 9z 




9U 3F 9R 9G 9R 3H 3R 

9t 9R 9r 9R 36 9R 9z 


Equation (14) is in weak conservation law form, while equation (17) is in 
non-conservation law form after the coordinate transforms have been intro- 
duced. The weak conservation law form may be retained, reference 5 and 6, 
as: 
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J Is the transform Jacobian 
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These equations may be rewritten in the following form in which the 
contravariant velocity components appear directly. 
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where 

(27) 

(28) 
(29) 

As described in reference 6, U , U and U are the contravariant 

R 9 X 

velocities perpendicular to the R, 0 and X coordinate surfaces. If these 
coordinate surfaces are chosen to coincide with physical boundaries in the 
problem, the no flow through solid surface boundary condition becomes; 
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Full conservation law forms (CLF) are generally to be preferred 
over non-conservation law forms (NCLF) on a theoretical basis. CLF contain 
the correct shock jump conditions while NCLF may not. The CLF, when used 
with appropriate spatial differencing, should have superior mass, energy and 
momentum conservation properties. Vftien a CLF finite difference solution is 
summed over a large volume of several mesh cells, the inter-cell flux of 
energy mass, and momentum cancel, thus satisfying the Integral conservation 
law properties for the large volume. For flow regions without shock waves, 
CLF solutions and NCLF solutions approach each other as the grid spacing 
decreases, but incorrect shock jump conditions may be predicted by NCLF solu- 
tions even in the limit of infinite grid resolution. 
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For simple flow problems where shock waves may be aligned with cell 
boundaries, the theoretical advantages of the CLF are easily demonstrated* 

For practical flow problems where shock waves may cross cell boundaries, the 
CLF advantages are more difficult to demonstrate. In fact, NCLF solutions 
often appear to give superior results. For this reason, the computer code 
described in this report has the capability to use either the NCLF of 
equation (17) or the weak conservation law form of equation (18). 
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Coordinate Transforms for Compressor Blade Rows 
As was previously discussed, the analysis will be limited to calcu- 
lations for the flow through a single blade passage. The physical flow 
domain is Illustrated in Figure 1 and consists of the space bounded by a pair 
of blades and the extension of their mean camber lines upstream and 
downstream of the blade row. The extension region is assumed to extend to 
upstream and downstream infinity, and the flow is assumed to be periodic, 
blade-to-blade , in this region. This physical domain is mapped onto a com- 
putational domain which is a rectangular parallelopiped . The computational 
domain is usually truncated 3 to 5 chords upstream and downstream of the 
blade row. 

Domain Regularizing Functions 

Two mapping functions were chosen: first, one which regularizes the 
physical domain; and second, one which locally increases mesh density near 
the blade leading and trailing edges. The domain regularizing functions are 
given by: 


R = 




(33) 
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Here is the radial position of the hub casing 

r^^p is the radial position of the tip casing 
0pg is the theta position of the blade pressure surface 
9 is the theta position of the blade suction surface 

S3 

z is the plan view projection of the blade row leading edge 
X* 3 

is the plan view projection of the blade row trailing edge 

A plan view for a typical compressor blade row in both physical 

space and computational space is shown in Figure 2. The hub and tip casings 

and the blade leading and trailing edge lines are mapped to straight lines . 

These mappings are adequate for nearly any compressor flow path. 

Similar views of the physical and computational space for blade to 

blade sections are shown in Figure 3. Here the curved blade shapes are 

mapped to untwisted planes. This mapping appears adequate for thin 

compressor blades but may not be adequate for thick blades with blunt 

trailing edges, as might be encountered in turbine blade rows. 

The effect of an axisymmetric part-span vibration damper or shroud 

can be included by modifying only the coordinate stretching in the r 

direction. Equation (33). A plan view of a typical flow path which includes 

a part-span shroud is shown in Figure 4. This geometry can be treated by 

defining two mapping functions; one function for the physical space below the 

damper (r < ^lower^ ’ ^ second function for the physical space above the 

damper (r > r ) . 

upper 
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^ lower ’^hub 


for r < r. 

- lower 


R_ = (1-C. ) u££er — ^ ^ 

2 1 r . - r 1 upper 

tip upper 


Here, ^ function of (r,z) specifying geometry of upper 

surface of damper, is a function of (r,z) specifying geometry of lower 

surface of damper, and is a constant determining placement of damper in 
computational space. The best value of .C^ is dependent on boundary con- 
dition formulations and will be discussed in a later section. This transform 
maps the damper into a slit in R-^ computational space as shown in Figure 5. 

If the blade leading and trailing edges are not constant z lines 
then near the upstream and downstream edges of the computational domain, the 
grid lines produced by these mappings are skewed with respect to either the 
local streamlines or a constant z line. This skewness may be eliminated by 
using the following transform, applied outside the blade row only, which 
allows a constant ^ line to approach a constant z line as z 4<». 
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*'®0 
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- 0.5 





19 


where ^teQ axial locations of the blade leading and trailing 

edge at the hub radius 

is a relaxation factor typically equal to 3.0 

2 * * z - upstream of the blade row, and 

z* » z^ “ z downstream of the blade row. 
te 

Mesh Packing Function 

A grid distributed uniformly in the (R,0,O space may be com- 
putationally inefficient since a fine mesh is often desired near the blade 
leading and trailing edges and a coarse mesh is desired near the upstream and 
downstream boundaries . A mesh packing function which satisfies both of the 
requirements is : 


£ = A sinh (K X) + B tan ^[sinh (K X) ] 

^ X X X X 


(39) 
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and is a free parameter which controls the mesh packing. Details on this 
transform are given by Merrington^. 
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This transformation from (R,0,5) space to (R,0,X) symmetrically 
packs points near the leading and trailing edges and stretches the domain 
away from the blade row. A typical grid structure having uniform distribu- 
tion in the (R,0,X) space is shown in the physical space in Figure 6 (plan 
view) and Figure 7 (blade-to-blade view) . As shown in these two figures , the 
grid system adopted is an offset one in which boundaries are located midway 
between grid lines. This arrangement has certain advantages in implementing 
both hard wall and periodic boundary conditions . These advantages will be 
discussed in a later section. These two figures also show that the blade-to- 
blade grid planes , R-^ or R-X planes , are not curved in the theta direction 
which means that the blade-to-blade lines are not normal to blade surfaces. 

In addition, dR/dd and dX/d0 are zero as may be seen from Equations (33) and 
(35). 
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Numerical Integration Scheme 
Equation (17) or (18) is solved using MacCormack's*" split operator 
finite difference scheme. This is a two step, explicit, second order 
accurate, conditionally stable scheme. The difference equations for the non- 
conservation law forms are: 
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u* = u“ - f 


7f" ,-f" 


M. ' 

[ j+l,k,£ j,k,£ 

3r n+l,k,i j.k.JlJ 

3z 


(44) 


U 


n+1 



' 

F* — F 

j,k,£ j— l,k,£^ 
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(45) 


where j, k, and £ respectively identify the axial, circumferential and radial 
coordinate of a grid point, n Indicates the time step, and * Indicates an 
Intermediate time step* Also 
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F* = F(U*) 

(46) 

H* = H(U*) 

(47) 

K* = K(U*) 

(48) 

= [O 1 0 0 0] 
R 

(49) 

Aq = [0 0 1 0 0] 

(50) 


Equations (40, 42 and 44) are the predictor steps, and equations 

(41, 43 and 43) are the corrector steps. The forward and backward space dlf-~ 

ferences appearing In the predictor and corrector steps are permuted each 

2 

operator execution as suggested by MacCormack. Similar operators are used 
for the conservation law forms. 

For this type of time splitting, the three MacCormack operators 

must be combined In S3rmmetrlc sequences In order to maintain second order 
2 8 

accuracy. * Denoting equations (41 and 42) as the Lj^(6t) operator, 

equations (43 and 44) as the LQ(6t) operator and equations (45 and 46) as the 
Ii^(6t) operator, a sequence may be written as: 

= [I^(5t) Lg(5t) I^(26t) Lg(6t) l^(5t)] u” (51) 

This simple but non-unique sequence advances the solution from time level n6t 
to time level (n+2)6t. 

In order to stabilize the solution along steep gradients of the 
flow quantities, such as across shock waves and near the leading and trailing 
edges , artificial viscosity terms were Introduced , that add to the numerical 
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damping of the scheme. The Lapldus^ form of these terms was used and for 
operator they are: 


Z. V 9 = 


i - u 


j,k,£+l j,k,£ 


V 


n n 

u - u 

''j,k,£+l ''j,k,i?. 

J ,k,£ J ,k,£-l J 


(52) 


where k is a parameter of order 1 that controls the strength of the artifi- 
cial viscosity terms. Low values of < produce sharp shocks, but result in 
unacceptable oscillations downstream of the shock. High values of < yield 
smooth solutions, but also very wide shocks and strong distortions to the 
inviscid field. A typical value of < is 0.4. 

Each one-dimensional operator is explicit and in the absence of 

artificial damping (Z^. , =0) have the stability conditions: 

J » k , 



(53) 


(54) 


( 55 ) 
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For calculations with strong shocks, these maximum time steps must be reduced 
by up to 40 per cent to maintain stability. 

For a given calculation 6t , 5t- and 6t_ may differ by over an order 
of magnitude, and the optimum value of 6t is selected using the principle 

that the 6t used with each operator be as close as possible to its maximum 

allowed value while maintaining a symmetric sequence. Typical values for 
transonic compressor blade rows are: 

= 20 (56) 

A suitable, symmetric sequence in this case would be 

(LQ(26t))\(20 6t)(Lg(26t))^ (57) 

where the notation (LQ(p6t))™ means m successive applications of the 

Lq operator with a time step of p times 6t. 
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BOUNDARY CONDITIONS AND SPECIAL TOPICS 

The computational domain Illustrated In Figure 8 Is bounded by open 
surfaces (flow may cross surfaces), solid walls and Inflow-outflow 
boundaries • Formulations and implementation of boundary conditions at each 
surface and Kutta condition are discussed In this section. 

Solid Wall Boundary Conditions 

The solid wall boundary condition appropriate for inviscid flow 
computations Is that of zero mass flux through the surface. This condition 
is difficult to Impose with uniform numerical accuracy in Euler equation 
simulations because of difficulties in evaluating derivatives of fluid pro- 
perties and velocities at the wall. In potential flow calculations this dif- 
ficulty does not appear since this boundary condition reduces simply to 
derivative of potential in direction normal to wall, being identically zero. 

A concise illustration of the problem is given by considering the 
integration of the continuity equation in two dimensions, as illustrated in 
Figure 9. The continuity equation is: 

If + ^ ° <58) 


or since v = 0 at wall 


+ p 


3v 

9y. 


9p , 3 N 

3t dx 


0 


0 


(59) 
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If the continuity equation Is to be numerically Integrated using 

3 I 

an explicit technique, the term ■=— (pu) provides no difficulties but the 
- - *0 

is difficult to evaluate with more than first order accuracy. 


term p 


9v 


3y 


0 


Using Taylor series: 


^i,l ^i,0 [3yJ 


5y + 


^^2 1 
9 V 


9y^J 


+ e?(6y^) 


(60) 


9v 


[9y 




f 2 ^ 
6 V 


6y 


6y^ 




( 61 ) 


No information about the second derivative of v is available to 

enable evaluation to second order accuracy. 

It is important to remember that this difficulty occurs for both 

the conservation and non-conservation law forms of the continuity and energy 

equations. In more complex flow geometries, where some type of body fitted 

coordinate system like that introduced in the previous section is used, this 

problem is not solved, only disguised. 

Several numerical flow solutions have been published using this 

reflection formulation with apparent good results (see references 10 and 11). 
12 

In addition Kreiss has demonstrated that for linear hyperbolic equations it 

th tti 

is sufficient for global n order accuracy to maintain n-1 order accuracy 
at the boundaries. The relevance of this statement for the Euler equations 
which are hyperbolic but nonlinear has not been established, but published 
solutions indicate that it may be at least a good guide for boundary con- 


ditions formulation 
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For simple flows or two-dimensional supersonic flows, the method 

13 

of characteristics may be used to compute flow solutions. Abbett has 
published a comparison of several boundary condition approximations coupled 
with second order accurate solution schemes, MacCorinack' s method in 
particular, witii method of characteristic solutions. Abbett found that 

usually provided 


dy 


0 


9v 

3y 


with 


0 


indeed the simple first order accurate scheme for 
more accurate answers than more complicated schemes to calculate 
second order accuracy. Another important result was that second order 
accurate approximations using one-sided difference expressions often gave 
very poor results. 

It is important to recognize that these problems with specifying 
the solid wall boundary conditions exist only for Integrating the fluid 
equations on the solid boundary. To integrate the fluid equations for in- 
terior points only the pressure on the solid boundary must be known. Solid 
boundaries, hub and tip casings or blade surfaces are mapped to coincide 
with coordinate surfaces which means that the contravariant velocity com- 
ponents U , U^, U (equations 30, 31, and 32) will be identically zero on 
these surfaces. MacCormack’s method or any scheme using centered spatial 
differencing will require either F or G, equation 24 or 26, on the boundary 
but not both. If = 0, F is determined only by the wall pressure and known 
coordinate derivatives. If = 0, G is also determined by the wall pressure 
and known coordinate derivatives. 

Fortunately the wall pressure can be calculated using the momentum 
equation in the direction normal to the solid surface. Following reference 
10, the normal momentum equation may be expressed as: 
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( 62 ) 


j »k,0 


where n is the surface normal, is the blade row rotation speed, and 
u are streamwise velocity components in the orthogonal boundary layer like 

T 

coordinate shown in Figure 10, and R^ are surface radii of curvature in 
this coordinate system, and (i^, i^) are unit vectors of the base (r,0,z) 
coordinate system- 

To be consistent with MacCormack's second order accurate scheme. 
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(64) 


(65) 


The term (pu /R ) can be approximated by 
s s 
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2 

and so be evaluated with accuracy ) if (3p/3n). and 

3 » ^ »1 

( 3u /3n) . are evaluated with accuracy i3'^(6n). These derivatives can be 

s 3 » 

evaluated with this accuracy by a number of schemes; simple extrapolation, or 

the first step of MacCormack' s method are two examples. The remaining two 
2 '' 

terms p(u^/R^) and ^^^ 2 ^^rel^ expressed similarly. 

An important aside at this point is that there is some theoretical 

basis, Kreiss reference 12, for the assumption that (3P/3n) need be evaluated 

with only accuracy c9^(5n). This evaluation, using (66) as a guide, requires 

only that we use (p) , f. = (p). ^ , (u ) = (u ) and (u ) 

J,K,U J.K.I, s j,k,o ® j.k.l j.k.o 

(u ) 

A comparative study conducted with reference 14 has shown this par- 
ticular scheme, with first order accurate (3P/3n) evaluation, to be the best 
available scheme for mixed supersonic and subsonic flows . For purely super- 
sonic two-dimensional flows, the method of characteristic type boundary for- 
mulations proved best. For purely subsonic irrotational flows, all schemes 
tested were adequate. The scheme proposed here performed well in both flow 
regimes, and for these two-dimensional flow tests almost never required arti- 
ficial damping for stability. 

While (3P/3n) is easily calculated numerically, the determination 
of the wall pressure requires that this pressure gradient be expressed in the 
stretched coordinate system (R,0,X). The normal pressure gradient is 
simply the dot product of gradient P and the unit normal vector. 

'3JP' 

3n 

V. J 


= (VP) • n 


(67) 
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The components of n, in physical space, are [y , Y-, y_] • For 

IT w Z 

simplicity the pressure calculation on the hub or tip casing will be 
discussed, which means that Yq= 0- The expanded form of equation (67) is 


^3nj 


[8rJ 







0, Y ) 
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or 


8n 


* Y. 


'apl 8R 
^8Rj 8r 
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3X 3r 


+ Y. 
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-T 

3X 

3Z 


C68) 


(69) 


Finally the pressure gradient may be expressed as: 




r 3r^ 


+ 



3R 

3Z 

J 




3X 

3r 


3X 
^Z 3Z 


IP 

3X 


(70) 


This expression for (3P/9n) may be considered as defining the components of n 
in the computational space. 


n 


comp 


3R ^ 3R 

3r 3Z 


. 0 , 


3X 


3X 


^ ^ ^Z 3Z 


(71) 


Equation (70) expresses the pressure gradient along one grid line, 
(3P/3R), in terms of the normal pressure gradient, (3P/3n), and the pressure 
gradient along a second grid line, (3P/3X). The pressure gradient (3P/3R) 
may be used to find the wall pressure. 
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Using the notation of Figure 9, the wall pressure, P, , » may be found as: 

1 > K , u , 


p. 

3 , 


k,0 


j ,k,l 


fdp' 

6r if 

f3p' 

[srJ 

j »k, 0 

9nJ 


is known to 


^(6n) 


(72) 


or 


■j,k,0 3 


4 P. 


- P. 


j,k,l jjk-,2 


’ 3 'f 


if 


is (73) 


or 


known to (^(6n ). Similar relations for the blade surface pressure may be 
derived . 


Problem Areas within the Standard Boundary Formulation 
The standard boundary formulation of equations (62) , (70) and (72) 
or (73) while complete suffers from several practical problems which have 
limited its full application. The difficulties are that all boundary points 
are now "linked" and must be solved for simultaneously, that pressure is ill- 
defined at corner points, and that surface radii of curvature, and R^ , are 
difficult to calculate accurately. 

The appearance of (3P/9X). , ^ in equation (70) is responsible for 
linking the boundary points together since it must also be approximated by a 
finite difference expression: 


9Xj^ 


'j,k,0 


^j+l,k,0 ^j-l,k,0 



2 A X 


(74) 


for example. This linking problem is quite vexing for a purely explicit cal- 
culation and has been dealt with by the expedient of lagging the calculation 
one time step so that (3P/3X) , „ at time level n is approximated by 

J 

(3P/3X). , - at time level (n-1). This approximation is quite good for the 

J > k» U 
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steady state or quasi-steady state problems generally encountered, and has 
also been used with good results in references 10 and 15. 

Corner points in a three-dimensional calculation present a dif- 
ficult problem in these calculations since multiple definition of pressure at 
these points can occur. The standard formulation provides little guidance on 
the proper choice since the bounding surfaces are not orthogonal. The impor- 
tance of proper definitions is illustrated by the disparity in reported solu- 
tions for external flow corner problems in references 16 and 17. The corner 
point problem has, in spite of its importance, been superceded by the problem 

of accurate calculation of surface radii of curvature R and R . An 

s T 

illustration of this problem is shown in Figure 11, which shows a typical 
supersonic blade section and the surface radius of curvature calculated from 
manufacturing coordinates using both spline fit procedures and finite dif- 
ference procedures. Curvatures are seen to be quite irregular and even 
discontinuous at x/C = 0.7. These specific problem areas plus general dif- 
ficulty in calculating second derivatives from manufacturing coordinates 
implies that the surface radii of curvature are known to first order accuracy 
at best. 

Since the curvatures are known, with present coordinate 

specifications, to only first order, the first order accurate equations for 

(dP/dn) and wall pressure are sufficient. The radii of curvature variations 

do of course cause oscillations in the wall pressure, but more importantly 

oscillations in (3P/3X), , » which is an important term in equation (70), In 

3 

blade tip sections these errors usually produce unacceptable variations in 
wall pressure and in extreme cases instability. When accurate curvatures are 
available (for example in the calculations presented in references 10 and 18, 
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these problems do not exist; but In order to proceed with calculations for 
less accurately known blade shapes, some further assumption Is needed. 

The best practical assumption appears to be evaluating (3F/3X) at 
location i. ^ 1 and not Jl * 0. This assumption, coupled with first order 
accurate evaluation of (9P/9n), Is equivalent to evaluating (9P/9n) on grid 
line 1-1 with the assumption that the streamline radii of curvature are 
equal to R and R . 

ST 

Boundary Conditions in a Staggered Grid System 

The flux preserving boundary formulations discussed above have been 
used with the weak conservation law form, equations (18), but a more natural 
grid system for the non-conservation law forms, equations (17) was suggested 
by Roache.^^ This system retains most of the advantages of the flux pre- 
serving boundary conditions while making the Imposition of periodic boundary 
conditions of the blade row and the Kutta condition much easier. This mesh 
system and the equations (17) are included In the associated computer codes 
as well as the weak conservation law forms and consistently give good 
results. This mesh system Illustrated In Figure 12 has no grid lines on the 
solid boundary, but maintains the solid boundary midway between an interior 
grid line and a dummy grid line . Meaningful fluid quantities are calculated 
only for points Interior to the flow and the exterior points are defined only 
to Implement the flux preserving boundary conditions. Wall quantities are 
defined as having average values of 


f . 

x,w 



(75) 


For the two dimensional example Illustrated In Figure 12, dummy values of 
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density, velocities and energy are found from reflection: 
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(76) 


The pressure is found as in the first mesh systems from equation (70). 
Following Roache,^^ these relations may be seen to be nearly equivalent to 
the flux preserving conditions. 

First. V = "T (v. o V. ,) *= 0 as desired. Then wall mass flow 
’ w 2 1,2 1,1 


(pv) = 
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p,- 


i>l i.l 
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1,1 w 


(77) 


Next wall momentum flux 
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(78) 


Thus this formulation does not directly conserve normal momentum flux. As 
Roache^^ points out, this flux term is still consistent with the original 
equation in the limit of (6x,6y) 0, Evaluating the momentum flux terra as 

would occur in any centered explicit scheme yields: 
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as (6y) -► 0, 


2 1 
Pv w+1-^ 

6y 


as desired, so that system is consistent. 


Periodic Boundaries 

One advantage of the staggered grid system is the logical simplifi 
cation it brings to periodic boundaries as well as to leading and trailing 
edge conditions . This grid system is illustrated in Figure 13 for the blade 
to-blade plane. Since no grid points are on the blade boundary, there are 
no stagnation streamline like grid lines which divide and go around a blade. 
Since the X and R grid stretchings are independent of 0 , the blade-to-blade 
running grid lines remain at constant (r,z) positions. The periodicity con- 
ditions are implemented simply by using the fact that the conditions at k=l 
are those of k=NTH-l and conditions at k=NTH are those of k=2 as shown in 
Figure 13. 

For the nonstaggered grid illustrated in Figure 9, periodicity may 
be maintained by: 
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If 




= 

j,NTH,A* 


the periodicity condition, then 


i *^j,NTH,Jr 


Upstream or Inflow Boundary Conditions 
Work described in this report has concentrated on isolated tran- 
sonic compressor rotors for which it is difficult to establish proper inflow 
upstream boundary conditions. It is well known that the bow shock system of 
such rotors produces an upstream disturbance that, while small, propagates 
far upstream as an acoustic disturbance. Two different studies have con- 
centrated on this problem, reference 19 for small disturbance theories and 
reference 20 for full Euler equation simulations. These studies adopted a 
common model, that of an isolated rotor operating in an infinite duct. The 
flow upstream is assumed to obey linear potential equations with the near 
field and the far-field potentials matched on a surface several chords 
upstream of the rotor. The upstream potential is expressed as the sum of 
three components: a uniform axial flow, a two dimensional axisymmetric per- 

turbation and a three dimensional perturbation. 


(p =(p + 0 . (r,Z) + (j) 

up uniform axi 


3D 


(r,0,Z) 


(82) 


Extensive simulations for one rotor geometry in reference 19 using non- 
reflecting boundary conditions showed that if the matching plane was located 
three or more axial chords upstream then setting , 0 ,Z) =0 was an 

adequate approximation for the near-field solution. These same simulations 
also demonstrated that the axisymmetric disturbance remained important even 
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three chords upstream. Nearly the same conclusions were reached in reference 
20 for two different configurations. 

These potential flow results suggest that if computational domains 
extend at least 4 to 5 chords upstream of an Isolated rotor then it is not 
necessary to model either the three-dimensional or axisymmetric disturbances. 
Computational domains of this extent greatly increase computation times and 
domains extending no more than 2 to 3 chords upstream are desired for reason- 
able computation times . Computational domains of this extent require that 
the axisymmetric disturbances be modeled with reasonable accuracy. 

A simple upstream flow model which can be used to simulate either 
potential or axisymmetric rotational upstream flows was proposed in reference 
10. This model represents the upstream flow as an unsteady axisymmetric base 
flow plus a more general but small disturbance. For the boundary conditions 
adopted for this report, the small disturbance will be neglected; in the 
potential flow case this is equivalent to neglecting <jj^j^(r , 0 ,Z) . The 
unsteady base flow is assumed to evolve from a steady axisymmetric 
"far-upstream" base flow, which is partially specified and partially computed 
from the interior solution. These boundary conditions are derived by repre- 
senting the velocity as: 

V = V(t,r,z) (83) 


the pressure as: 


P = P(t,r,z) 


(84) 
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and the density as: 


P = p(t,r,z) 


(85) 


Note that this definition will include both 
19 in the base flow V and P . 


^uniform 


and 


d) . from reference 
axi 


o o 

The equations describing the base flow are continuity and the three 
momentum equations and continuity; 


and 


DV _ _ Vp 
Dt " P 


§1 + (PV) = 0 


( 86 ) 

(87) 


The base flow is also assumed to be isentropic. 

Using as a coordinate system the modified cylindrical coordinate 
system shown below, the equation of motion may be written in terms of the 
tangential velocity and the meridional velocity, 





The equations of motion written in a general orthogonal curvilinear coordinate 
system are 
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Since = 0, by definition and = 0 by assumption, these 
equations reduce for the base flow to, 
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The first two equations may be rewritten in characteristic 

o 

form using the isentropic relation, dp = a dp , as 
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Adding and subtracting these equations yields 
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The third equation, (93), is already in characteristic form. 
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A fourth characteristic also exists and has been chosen to be 
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5l=° 


( 101 ) 


which Implies that an Isentroplc flow exists along each streamline from the 
”far-upstream" condition to the computational domain inflow boundary. 

At each point on the inflow boundary four pieces of information 
must be supplied in order to calculate the flow variables. For transonic 
compressor rotors, the upstream meridional Mach number is always less than 
one, and the J characteristic is an upstream running characteristic. The 
J characteristic is always a downstream running characteristic. The charac- 
teristics associated with equations (100) and (101) are simply the convection 
path. Thus three pieces of information may be specified independently at the 
inflow boundary. 

This situation is Illustrated for a single spatial dimension in 
Figure 14. For this formulation the solution does not exist for x < 0 and is 
assumed to be known at time level n6t. At the node point (n+1, 1) three 
characteristics are incident: a convection characteristic, an upstream 

running characteristic, and a downstream running characteristic. The 
value of the characteristic is known, but the information carried along 
the remaining characteristics nuist be specified. 

r+ ^n+l [ 'in+1 , 2 n+1 


(J 


t- .n-t-j. f 'ini-x , z nH 
r \ ^ -1 


( 102 ) 


(J‘ 
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At the open boundary (j=*l), (u^)*^^^ or may be specified in order 

to determine the characteristic value. When is specified, 


takes the value: 
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+ t“ . ^ 

J = J . + — ^ a_ 
r 1 y-1 1 


When (u )“ is specified, J takes the value 
nij^ r 


J-^ = 2(u J- 

r m ^ 1 


(104) 


(105) 


For the computer codes described in the report, the inflow boundary 
is assumed to be located far enough upstream that the radius of curvature of 
streamlines is large enough that the curvature term (l/h^h^ )(8hQ /3m) may be 
neglected. The resulting equations describing the upstream flow are: 


D J 
Dt 

= 0 

(106) 

d'j” 

Dt 

= 0 

(107) 

CD 

= 0 

(108) 

Ds 

Dt 

= 0 

(109) 


The constraint equation, equation (94), is not used and the grid is 
constructed such that the hub and tip casing slope are zero at the inflow 
boundary so that u^ may be assumed to be zero. 

Two input options are provided in the computer code for specifying 
the inflow boundary condition. The primary option uses the freedom to spe- 
cify a^ to maintain constant values of total pressure and total tem- 
perature at the inflow boundary. The user may specify the values of total 
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pressure and temperature as a function of radius but may not specify the 

Inlet velocity. This option Is the normal one and does not specify the Inlet 

mass flow rate. With some Initial conditions, this Inflow boundary condition 

option can become unstable. A second option Is provided for use with the 

unstable option 1 cases. In option 2, the value of the downstream ninnlng 

characteristic, J^, Is specified directly by specifying the desired value of 

the uniform far-upstream Mach number. The value of the characteristic is 

r 

calculated as: 



r 


(n ) 
m up 



a 

up 




( 110 ) 


Use of option 2 will usually lead to a small total pressure error at the inflow 
boundary. 

For both options, once values for and are known, the values 
of and are calculated using equations (102) and (103). The 

values of T^ , ^ calculated using the Input values of 

total temperature and the Isentropic flow relations. The value of 
(Ug)^^^ is specified by the user and the radial velocity (u^)^^^ Is 
assumed to be zero. For option 2, the characteristic value is specified 
by the user, while, for option 1 the characteristic value is selected 
through iteration on the value of a“"^^ to match the input value of total 
pressure. 

Downstream or Outflow Boundary Conditions 
Outflow boundary condition formulations closely follow the model 
for the inflow boundary conditions, equations (91) through (94). However, 
for subsonic flow, as shown in Figure 15, only one characteristic, j“, 
brings Information Into the computation domain. For supersonic flow, no 
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characteristics bring Information Into the computational domain. For sub- 
sonic flow, the information to be specified is the far-downstream static 
pressure, while, for supersonic flow, no downstream boundary condition may be 
Imposed. The downstream boundary Is assumed to be located far enough 
downstream that the streamline radius of curvature may be assumed to be zero 
and the radial velocity may be assumed to be zero. The equations to be 
solved are: 


and 


) J 
Dt 

= 0 

(111) 

Duq 

Dt 

= 0 

2 

(112) 

3P 


(113) 

3r 

r 


Dt 

= 0 

(114) 


Since the swirl velocity on the outflow boundary is determined by the 
interior solution, the exit static pressure may not be specified arbitrarily 
as a function of radius. To avoid this problem a static pressure on the 
inner casing is specified and the radial static pressure variation is calcu- 
lated from equations (113) and (114) . The exit static pressure variation is 

calculated by first extrapolating the theta-averaged values of and u„ to 

X 

the last computational plane (j“Nx) from the upstream computational planes 
( j=NX-l and j=Nx-2) . The resulting centrifugal force gradient at the outflow 
plane is balanced by the radial pressure gradient. The radial pressure gra- 
dient is calculated from equation (113) , using the extrapolated values of 

u- and density, p, by trapezoidal rule integration. The final exit density 

y 

and speed of sound, a, are then calculated using the isentropic flow 
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assumption, equation (114). The exit velocity is then calculated from the 
value of the extrapolated characteristic. The radial velocity is also 
assumed to be zero. 

These boundary conditions are equivalent to the physical 
constraints imposed when a compressor is tested. Inlet or reservoir con- 
ditions are specified, usually atmospheric stagnation conditions, and the 
inlet swirl velocity is specified, usually zero. Downstream of the 
compressor, some throttling device is used to determine the compressor 
operating point. A particular operating point is determined by the intersec- 
tion of the compressor characteristic and the throttle characteristic. The 
operation of this device could be simulated by specifying either the 
compressor static pressure ratio or the mass flow. It is computationally 
somewhat easier to specify the rotor static pressure ratio. 

No attempt was made to calculate the actual flow at the throttle or 
exit boundary, and hence account for the effect of disturbances there. Such 
a calculation, especially in the strongly rotating flow downstream of the 
rotor, is potentially as difficult as the through-flow calculation itself. 
This simplification does imply some inconsistencies at the boundaries, while 
at the upstream boundary the transmission of acoustic disturbances is 
falsified; at the downstream boundary vortlcity convection is distorted. The 
acoustic disturbances carry very little energy and their falsification has 
been shown to have very little influence on the flow over the rotor. When 
the downstream perturbations are not calculated, the effective model of the 
flow beyond the downstream computational boundary is one which is axisym- 
metric and uniform in the streamwise direction. Such a flow can have no 
radial vortlcity and is potentially inconsistent with the blade flow solution 
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if that solution requires a significant amount of radial vorticity to be pre- 
sent at the computational boundary. For the present calculations, this 
inconsistency does not appear to be important even though some radial vor- 
ticity is shed by the rotor. 

Kutta Condition 

The Kutta condition, a uniqueness condition Imposed on Inviscid 
flow solutions in order to approximate the true viscous flow solution, is 
imposed quite simply for the compressor rotors studied. No attempt is made 
to model blunt or cusped trailing edges. Instead, the flow at the last axial 
grid points on the pressure and suction surfaces is assumed to be parallel to 
the blades and to have equal static pressures. This condition is adequate 
for compressor rotors designed to have subsonic outflow, but would require 
modification for turbine blade rows having supersonic outflow. 

Spinners 

Several attempts were made to model a spinner, including one that 
would smoothly merge with the centerline, but these attempts were 
unsuccessful. High radial velocities tended to build up at the points close 
to the centerline, resulting in mass flow defects of the order of 5%. Since 
r=0 is a singular point in a cylindrical coordinate system, the implemen- 
tation of boundary conditions at the centerline can pose some problems . In 
addition, the centerline is a stagnation streamline, and stagnation points 
pose stability problems. In view of the above problems, the spinner was 
finally replaced by a cylindrical centerbody with a hub to tip ratio of 0.3. 
Several promising attempts to deal with centerline procedures in a cylindri- 
cal coordinate system have been published (see reference 21 for example), and 
the spinner problem should now be re-examined. 
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COMPUTER CODE DESCRIPTION 

The numerical procedures previously described have been implemented 
in two computer programs. The first, MESH3D, generates the finite difference 
grids, interpolates blade geometry data from manufacturing coordinate sec- 
tions to the finite difference grid, computes blade normal vectors and 
curvatures, and generates coordinate stretching derivatives (3R/9r) etc. 

The second, BLADE3D, generates initial conditions if desired, integrates the 

appropriate equations of motion, and prints a solution matrix. A third pro- 
gram, GRAPH3D, produces two output files which may be used for user defined 
graphics ^nd produces some rudimentary graphics output. 

Grid Generation Program: MESH3D 

The grid generation program has four important geometric inputs: 

1) (r,z) coordinate pairs describing hub and tip casing shapes; 2) (r,z) 
coordinate pairs describing location of blade row leading and trailing edges; 
3) (r,z) coordinate pairs describing damper geometry where applicable; and 4) 
(r,9,z) coordinates describing the blade row geometry. 

Operation of MESH3D on these geometric quantities is best described 
using Figure 16 (the axial grid numbering scheme) and Figure 17 (simplified 
MESH3D flow chart). The input to the main program is the grid numbering 
information, particularly the axial grid numbering scheme, which is necessary 
to size the finite difference grid and control its generation. Subroutine 
XRGRID controls the input of all axisymmetric geometric data and generates 
(r,z) coordinates of the finite difference grid and calculates some of the 
coordinate stretching derivatives [(3X/3z), (3X/3r), (3R/3r), (3R/3z)]. 
Subroutine GALT controls input of blade surface geometry, interpolates this 
data onto the finite difference grid, calculates blade normal vectors, and 
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calculates blade radii of curvature. Subroutine BLDT calculates [(30/3r), 
(30/30), and (30/3z)] and creates an input file for BLADE3D containing all 
the required blade geometry information. Subroutine BUILD calculates hub 
and tip normal vectors as well as the radius of curvature along the axial 
direction, and subroutine SVGRID creates input files for BLADE3D containing 
all required information about grid coordinates • 

Numerical Integration Program: BLADE3D 

The numerical integration program is modular in operation, can be 
executed either on large mainframe computers or on minicomputers with a mini- 
mum of code changes, and is inherently restartable. A simplified flow chart 
of BLADE3D is shown in Figure 18. The main program input describes the 
axial grid numbering scheme (which is identical to that of MESH3D) , describes 
the operator sequence to be executed , the integration time step , and the 
number of operator sequences to be executed. Subroutine OPEN begins the 
calculation setup process by reading two data files generated by MESH3D 
which describe the blade row geometry. Subroutine START completes the setup 
process either by reading a previous solution matrix (solution vector U at 
all grid points) or generating a new starting solution. Subroutine THREE 
controls printing of starting and final solutions as well as execution of the 
requested operator sequence. The finite difference operators equations 
(40-44) are implemented in subroutines BLOCKX, BLOCKT, and BLOCKR, 

Subroutine CLOSE outputs the final solution matrix. A lengthy calculation 
will usually be run in several sections , and this operation of START and 
CLOSE provides a needed restart capability since the intermediate results are 


available for backup 
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Graphical Output Program: GRAPH3D 

The graphics output program produces two solution matrix data 
files which contain all information needed for user defined graphics 
packages. The first file contains for each node point in the finite 
difference grid its physical space r,9,z coordinates, the solution vector 
[rp, rpu^, rpUg, ^ modified solution vector > 

T^.]* The second file contains the computed streamline positions inside 
the blade row volume. The starting point for each streamline calculation is 
a finite difference grid node at the rotor leading edge. This second file 
contains the (r,G,z) coordinates of a streamline position and the solution 
vector at that point, [rp, rpu^, rpug, rpu^, rE^, P^, M^, Mq, T^]- 

The computed streamline positions thus correspond to the traditional SI 
and S2 streamsurface definitions. 

In addition to these data files, GRAPH3D also produces printer/ 
plotter plots of blade surface pressure and Mach number for each radial 
grid plane and each S2 surface. 

For user convenience, a subroutine, called USERREAD, is included 
with GRAPH3D which may be used to read the two solution matrix data files. 
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INPUT DICTIONARY FOR MESH3D 

Input variables may be described in any consistent set of units 
since the MESH3D and BLADE3D programs contain only non-dimensional units. 
The length scale selected is usually the tip casing radius at the farthest 
point upstream; the velocity scale selected is the stagnation temperature 
speed of sound at this same point. Each time geometric quantities are 
input, a scaling factor is required which converts any dimensional 
geometric quantity into a nondimens lonal quantity. In this way, the most 
convenient set of units may be input for each geometric quantity, which is 
then nondimensionalized by the code. In addition when axial coordinates 
are input, a coordinate system origin adjustment factor is also required. 
This additive factor corrects possible coordinate origin differences bet- 
ween sets of geometric input . The input variables for MESH3D in the 
order they appear in Table 4.1, are the following: 

Main Program Input 

TITLE Title for problem identification, any information may appear in the 
first 70 columns. 

NBL Number of blades in row, must be 

NX Number of axial grid planes, must be 3 ^ NX £ 100 

NTH Number of theta grid planes, must be 3 ^ NTH ^ 17 

NR Number of radial grid planes' must be 3 _< NR _< 18 

*N1 Axial plane containing the blade row leading edge 

*N2 Axial plane containing the blade row trailing edge 

*Note: N2-N1 must be an even integer ^50. 
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KUTTA 

ILE 

ITE 

SPHALF 

IAN 

IWR 

lEMP 

KX 

A1 

TOL 


Axial plane at which Kutta condition is to be applied. 

KUTTA = N2 normally. 

Axial plane containing the damper leading edge. A value of 
0 indicates no damper. 

Axial plane containing the damper trailing edge. A value 
of 0 Indicates no damper. 

Mesh type parameter. If = 0, then first mesh system is 
generated. If =0.5 the second or staggered mesh system is 
generated. Only the staggered mesh may be used with the 
current version of BLADE3D, so SPHALF =0.5 should be used. 

Type of transform derivative calculation. If IAN = 0 
then metric calculation done by finite difference method 
using node (r,0,z) positions and if = 1 metric calcula- 
tion is done by analytic functions. IAN == 1 is 
recommended. 

FORTRAN unit number for MESH3D bulk output, may be equal 
to 6 or greater than 9. 

MESH3D debug output flag, allowable values are 

0 => none (recommended value) 

1 => SUBROUTINE GRID output 

2 => SUBROUTINE INTER output 

3 => SUBROUTINE GETTHA output 

4 => SUBROUTINE CALDR output 
5,6 => SUBROUTINE CALDER output 

-1 => All subroutines (not recommended) 

Axial coordinate packing factor. KX =0.1 corresponds 
to non-packed grids. KX = 3.0 corresponds to highly packed 
grids. See special instruction number 1 for aid in 
selecting KX. 

Radial Grid line relaxation factor, used to make com- 
putational space grid lines approximate Z = constant lines 
near inflow and outflow boundary. A value of 3.0 appears 
to be best. See special instruction number 1 for aid in 
selecting Al. 

Convergence tolerance on axial grid plane positions. A 
value of O.OOOl is a tight tolerance and a value of 0.001 
is recommended. 
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lAB Convergence choice parameter* A value of 0 will cause 

axial plane position calculation to abort after 10 
iterations. Nonconvergence is usually caused by input 
errors in either tip and hub casing positions or leading 
and trailing edge positions. A value of 1 allows uncon- 
verged plane positions to be accepted. lAB = 0 is 
recoramenided. 


Subroutine GRID Input 


HUB Coordinates 

NRH Number of coordinate pairs 

2 j< NRH < 50. 

SCALEl Input hub coordinates will 

(XH + ADJl) /SCALEl 
and 

RH = RH/SCALEl 


ADJl 


describing hub casing contour, 
be scaled as: XH = 


XH and RH Coordinate Pairs describing hub casing contour. 

TIP Coordinates 

NRT Number of coordinate pairs describing tip casing contours. 

2 ;< NRT _< 50. 

SCALE2 Input tip coordinates will be scaled as: 

XT = (XT + ADJ2)/SCALE2 
and 

RT = RT/SCALE2 


ADJ2 • 

XT and RT 

Leading Edge 

NZLE 

SCALES 

ADJ3 

ZLE and RLE 


Coordinate pairs describing tip casing contour. 
Coordinates - See Also Special Instruction Number 2 


Number of coordinate pairs describing leading blade edge 
location. 2 ^ NZLE ^ 20 

Input leading edge variables scaled as: 

ZLE = (ZLE + ADJ3) /SCALES 
and 

RLE = RLE/SCALE3 


Coordinate pairs describing blade leading edge contour. 
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Trailing Edge Coordinates 

NZTE Number of coordinate pairs describing blade trailing edge 

contour. 2 < NZTE < 20 

SCALE4 Input trailing edge variables scaled as: 

ZTE = (ZTE + ADJ4)/SCALE4 


and 


RTE = RTE/SCALE4 

ADJ4 

DAMPER COORDINATES - SKIP IF BLADE ROW HAS NO DAMPER 


Damper Lower Surface 


XTE, ZTE 
NZL 

SCALE5 

ADJ5 
ZLL; RL 


Coordinate pairs describing blade trailing edge contour- 

Number of coordinate pairs describing damper lower 
surface. 2 ^ NZL ^ 20 

Input damper lower surface variables will be scaled as: 
ZLL = (ZLL + ADJ5) /SCALES 
and 

RL = RL/SCALE5 


Coordinate pairs describing damper lower surface. 


Damper Upper Surface 


NZU 

SCALES 


ADJ6 
ZU, RU 
Cl 


TITLE2 

NPTS 


Number of coordinate pairs describing damper upper surface. 
2 _< NZU _< 20 

Input variables will be scaled as: 

ZU = (ZU + ADJ6) /SCALES 
and 

RU = RU/SCALES 


Coordinate pairs describing damper upper surface contour. 

Location of damper in computational space. 0 < Cl < 1.0; 
Cl may take only certain values, see special instruction 
number 3. 


Subroutine CALT Input 

Title Card for blade geometry input set. 

Number of points on a specification line, see Figure 19 
and special instructions. Maximum value is 25. 
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Number of stacked blade sections, see Figure 19. Maximum 
value Is 15. 

Input blade data scaled as: 

Z = (Z + ADJ7) /SCALE? 

R = R/ SCALE? 

Blade geometry Input flag, either 0 or 1. 

See discussion below for option forms. 

Number of blades in blade row. 

Coordinate origin adjustment factor. 

Three dimensional coordinate array of blade pressure and 
suction surface positions. First index varies from 1 to 
NL, while second index varies from 1 to NPTS. For the 
third index a value of 1 is used to store radial 
coordinate; a value of 2 is used to store axial coordinate; 
a value of 3 is used to store theta coordinate of pressure 
surface; a value of 4 Is used to store theta coordinate of 
suction surface. See also special instruction number 4. 

Variable 
MID(1,J,1) 

MID(NL,J,1) 

MID(1,J,2) 

MID(NL,J,2) 

MID(1,J,3) 

MID(NL,J,3) 

MID(1,J,4) 

MID(NL,J,4) 


Comments 

NPTS Z coordinate values for blade section #1 

NPTS Z coordinate values for blade section //NL 

NPTS R coordinate values for blade section //l 

NPTS R coordinate values for blade section //NL 


NPTS values of either or 9 for blade 

section //I 

NPTS values of either or 0 for blade 

section //NL ^ 


NPTS values of either 0 or 0 for blade 
section //I 


NPTS values of either 0 or 0 for blade 
section //NL 


NL 

SCALE? 

IM 

NBL 

ADJ? 

MID(I,J,K) 
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Special Notes and Restrictions on Input Variables for MESH3D 

1. The best values of KX and A1 can only be determined by trial 
and error. A small value of KX (0.01) leads to nearly unpacked X-R grids, 
see figure 20. Larger values of KX progressively pack X-R grid planes near 
the blade leading and trailing edges while increasing the grid spacing in the 
far field, see figures 20 and 21. The best value of A1 depends on KX and 
should be selected after KX is selected, see figure 22, 

2. Due to the placement of X-0 grid planes outside the normal hub 
and tip casings, it Is imperative that definitions of blade leading edge and 
trailing edge geometry be continued at least 10 percent of blade span outside 
the hub and tip casings. 

3. The constant determines the part-span shroud placement in 
computational space. The restriction on shroud placement is that it must lie 
midway between two X-0 grid planes. Thus the values of are restricted to 
FLOAT [( IDL-1) / (NR-2) ] . IDL is the number of the X-0 grid plane immediately 
below the damper. It is suggested that to generate a grid containing a part- 
span shroud that first a grid without a part“span shroud be generated. The 
values of ILE, ITE, and IDL may then be determined by inspection. 

4. Blade geometry is given by four two-dimensional arrays, each 
row of which describes a blade section. Blade sections are numbered from 1 
to NL. The first blade section at the hub and the last one at the tip shall 
not necessarily conform to the hub or tip casing profile, but may be given 
within the flow region, crossing the boundary, or completely outside of the 
boundary. Column 1 and column 2 are used to describe the axial and radial 
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positions of the blade section, while column 3 and column 4 will describe 
either the angular position of the mean camber line (in radians) and the 
tangential thickness (in radians) of the blade section or the angular posi- 
tion (in radians) of the blade pressure surface and the angular position of 
the suction surface of same blade* Option flag IM = 0 requires mean camber 
line and thickness while IM = 1 requires positions of pressure and suction 
surfaces. See Figure 23 for definition of tangential thickness. These blade 
coordinates are stored in the single three-dimensional array MID(I,J,K). 

5. All variables allocated 5 columns in Table 4.1 are 
integer variables and must be right justified. All variables allocated 10 
columns are real variables which should have decimal point included. 
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/ Hub Casing Coordinates 


XH(NRH) RH(NRH) 


XT(1) RT(1) 


XT(NRT) XT(NRT) 


NZLE SCALE 3 I ADJ3 


V Tip Casing Coordinates 


ZLE(l) 1 RLE(l) 



ZLE(NZLE) 


^ Leading Edge Coordinates 
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*** NOTE NEXT SET OF DESCRIPTOR CARDS *** 


t 


*** ARE PRESENT ONLY IF A DAMPER IS *** 


*** INCLUDED IN THE SOLUTION 
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Jfi 1.!^. t... 


MID(1,1/1) to MID(NPTS,1 ,1 ) NL SETS OF NPTS POINTS. BLADE Z COORDINATES. 




MID(1,2,1) to MID(NPTS,2,1) 





I 


MID(1,NL,1) to MID (NPTS, NL,1) 





5 MID (1,1, 2) to MID (NPTS, 1,2) NL SETS OF NPTS POINTS. BLADE R COORDINATES, 





MID(1,NL,2) to MID (NPTS, NL, 2) 



MID (1,1, 3) to MID (NPTS, 1,3) NL SETS OF NPTS POINTS. BLADE OR 

„ _ _ _ „ COORDINATES 

' ~ T ^ I I '] 'T" 

! II 


MID(1,NL,3) to MID(NPTS,NL,3) 


BLADE t or 0^^ 

MID(1,1,4) to MID(NPTS,1,4) IJL SETS OF NPTS POINTS. COORDINATES 


MID(1,NL,4) to MID(NPTS,NL,4) 
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INPUT DICTIONARY FOR BLADE3D 


The input variables for BLADE3D, in the order they appear in 

Table 4.2, are the following: 

TITLE Title for problem identification, any information may 

appear in the first 70 columns 

NBL Number of blades in row, must be 1 

NX Number of axial grid planes, must be 3 ^ NX ^ 100 

NTH Number of theta grid planes, must be 3 _< NTH ^ 17 

NR Number of radial grid planes, must be 3 ^ NR ^ 18 

N1 Axial plane containing the blade row leading edge 

N2 Axial plane containing the blade row trailing edge 

KUTTA Axial plane at which Kutta condition is to be applied. 

KUTTA - N2 normally. 

ILE Axial plane containing the damper leading edge. A value of 

0 indicates no damper. 

ITE Axial plane containing the damper trailing edge. A value 

of 0 indicates no damper. 

IDL Eladial grid line below damper, see Figure 16. If no damper 

then IDL = 1. 

IDH Radial grid line above damper, see Figure 16. If no damper 

then IDH = NR. 

IFCL IFCL=0 gives non-conservation law form equations, see 

equation 17, IFCF=1 gives full conservation law form, see 
equation 18, 

IROT If IROT=0 the fully unsteady equation sets, equation 17 or 18, 

are used. When IR0T=1 a pseudo-unsteady equation set is used 
in which the energy equation is replaced with a constant 
rothalpy assumption. 


GAMMA 


Ratio of gas specific heats, usually 1.4. 




61 


CAPPA 


PDOWN 


INOPT 


NUTH 


RUTH, UTHIN 
NTT 


RTT, TTIN 


NPT 


Blade row angular rotation speed parameter equal to 
(Wrj-ip/aQ). is rotation speed in radians per second, 
r is for upstream tip radius (reference length for non- 
dimensional lengths), a is far upstream stagnation speed 
of sound on tip streamline (reference velocity for non- 
dimensional velocities). 


Artificial viscosity parameter, usually equal to 0,5. A 
value of 0.1 would be very low damping while a value of 1.0 
would be high damping. 

The static pressure to be set on the last axial grid point 
at the hub radius. This pressure, (P/p a ^), controls the 
blade row pressure ratio. 

IN0PT=1 corresponds to inflow boundary condition option 1 in 
which the user specifies the upstream total pressure and 
temperature. IN0PT=2 corresponds to inflow boundary condition 
option 2 in which the user specifies upstream total tempera- 
ture and the values of the characteristic (see inlet B.C. 
discussion) . 

Number of radial positions and absolute tangential velocity 
pairs which specify inlet tangential velocity as a function 
of radius. If NUTH=1, a constant value of UTHIN is produced 
and the input value of RUTH is ignored. NUTH must be ^20, 

Radial position and absolute tangential velocity theta averaged. 

UTHIN = U^/a and RUTU is the scaled radius, RUTH = R/R 

0 o ref 

Number of radial positions and total temperature pairs which 
specify inlet absolute total temperature as a function of 
radius. If NTT=1, a constant value of TTIN is produced and 
the input value of RTT is ignored, NTT must be 5 20, 


Radial position and absolute total temperature, theta 
averaged. TTIN is the ratio between the local total tem- 
perature and the reference te mpera ture . The reference tem- 
perature is defined by a = / yRT and RTT is the scaled 

radius, RTT - ° ^ 

Number of radial position and total pressure pairs which 
specify inlet absolute total pressure as a function of 
radius. If NPT=1, constant value of PTIN(J^) is produced 
and the input value of RPT is ignored, NPT must be -20. 
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RPT, PTIN 
or 

RPT, JRIN 

UPMACH 

DT 

NSTEP 
IB EC 

NOR(I,J) 

lOP J = 1 

IML J = 2 
lEX J = 3 

NPl 

CODE(NPl) 

JST(NPl) 

JEND(NPl) 


Radial position and absolute total pressure, theta averaged. 
PTIN is the ratio between the local total pressure and the 
reference stagnation pressure. P = 1/y and RPT = R/R 

When INOPT *= 2, values are input^®^ instead of PTIN, 


Far upstream uniform absolute meridional Mach number. Used 
only when initial conditions are generated by BLADE3D. 

Scaled time step to be used. If DT = 0.0 then time step 
IS selected by program. See special instruction for rules 
for selecting DT, 


Number of split operator cycles to be run. Several hundred 
operator cycles are usually required for a solution. Each 
cycle advances the solution several time steps . Solution 
is usually run only a few hundred cycles at a time to pro- 
vide a restart capability. 

If IBEG = 1, a starting solution is generated. 


The MacCormack operator sequence to be run is stored in a 
two-dimensional array, NOR (10,3). Each row specifies an 
operator (X direction, 0 direction or R direction), the 
number of times each operation is to be executed and the 
time step multiple to be used. The entire operator 
sequence is executed NSTEP times. 


The individual codes associated with each row are; 

Operator Code 1+X operator, 2+0 operator and 3+R 
operator 

Time step multiple 
Execution time multiple 

See special instructions and input card format sheet for 
more information about operator sequence. 

Number of different sets of solution matrix information to 
be printed before operator sequence is begun. 

Determines type of printed solution variables. See special 
instructions for definition of allowed values. 

First axial station to be included in printed information. 

Last axial station to be Included in printed information. 
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LST(NPl) 

LEND(NPl) 

NP2 

NTIMES 

C0DE(NP2) 

NST(NP2) 

JEND(NP2) 

LST(NP2) 

LEND(NP2) 

NP 3 

C0DE(NP3) 

JST(NP3) 

JEND(NP3) 

LST(NP3) 

LEND(NP3) 


First radial station to be Included in printed information. 
Last radial station to be included in printed information. 


Number of different sets of solution matrix Information to 
be printed during operator sequence. If NP2 “0, no inter- 
mediate results are printed , 


Solution matrix is printed before the first operator 
sequence, after NSTEP sequences, and after every NTIMES 
sequence during the solution integration. NTIMES = 0 gives 
no intermediate printed results. 

Detemnines type of printed solution variables. See special 
instructions for definition of allowed values. 


First axial station to be included in printed information. 
Last axial station to be included in printed information. 
First radial station to be included in printed information. 
Last radial station to be included in printed information. 


Number of different sets of solution matrix information to 
be printed after operator sequence is completed. 

Determines type of printed solution variables. See special 
instructions for definition of allowed values. 

First axial station to be Included in printed information. 

Last axial station to be included in printed information. 

First radial station to be Included in printed information. 
Last radial station to be included in printed information. 
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Special Instructions for BLADE3D Input 

1. DT is the non-dimensional integration time step. DT = 

[ (5 ta^)/ ] . The maximum DT in the absence of artificial damping or 
shock waves is given by equations 53, 54 and 55 and is printed out with each 
execution of matrix print routine. This maximum value must be reduced by 
approximately 40% for calculations with strong shocks. The suggested maximum 
value of DT is printed at the end of each code run. 

2. The sequence of MacCormack split operators to be run must be 

determined for each case using the principle that (DT)^, ^^^^0 

should be as close to the maximum allowed value while maintaining a symmetric 

2 

sequence, as described by MacCormack. See also equation 51. 

Denoting the operators as L^(DT^), (DT^), then a simple 

sequence may be written as: 

which advances the solution from time level n to level n + 2. This sequence 

assumes (DT ) ~ (DT ) and (DT ) ~ 0.5(DT ) . A more typical 

X 0 X R 

max max max max 

sequence for transonic compressors arises when 

(DT ) ~ 10 (DTq) ~ 20 (DT^) 

iunx max max 

A suitable, but non-unique, sequence for this case would be: 

= [10 L^(DT^), 5Lg(2DTg),L^(20Dy, 5Lg(2DTg), 10L^(DyjU^ 



65 


In general an operator is described by three parameters: the base time step 
(DT), a time step imiltiple (1, 2 or 20 above), and an execution multiple (10, 
5 or 1 above)* 

The full operator sequence is stored in N0R(I,J) which for above 
example would appear as: 



Operator 


1 = 1 

1 = 2 

1 = 3 


Number 


Operator 

Time Step 

Execution 




Key 

Multiplier 

Multiple 

J 

= 1 


1 

1 

10 

J 

= 2 


2 

2 

5 

J 

= 3 


3 

20 

1 

J 

= 4 


2 

2 

5 

J 

= 5 


1 

1 

10 


The operator keys are l+L, 2->-L_, 3-»-L„ and the maximum number of opera- 

X o K 

tors in a sequence is 10. 

4. The value of CODE(I) determines the solution variable to be 
printed. The allowed values are: 
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Value Printed Variables 


1 Conservation variables (rp, rpu^, rpE^) 

2 Physical variables (p , u^, Uq , u^, E^) 

3 Physical variables in coordinate system rotating with 
blade row (p , M , M- , M , P) 

rel 

4 Physical variables in coordinate system rotating with 

blade row (P^, M , M_ » M > T^.) 

t* r* 0 ’ z* t 

rel 

5 Physical variables in laboratory coordinate system 

(p^. Mg, T^) 

In addition, when codes 4 and 5 are selected, a summary of mass-flow-weighted 
axisymmetric variables is presented. 

5, All variables allocated 5 columns in Table 4,2 are Integer 
variables and must be right justified. All variables allocated 10 columns 
are real variables and should have decimal point included. 
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Table 4.2. Input Card Format for BLADE3D 


1 

— 

6 

11 

16 

21 

TITLE 

NBL 





NX 

NTH 

NR 



N1 

N2 

KUTTA 

ILE 

ITE 

IDL 

IDH 

IFCL 

IROT 



26 


31 


36 


41 


46 


51 


56 
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GAMMA 


INOPT 


NUTH 


W 


CAPPA 


PDOWN 


RUTH(l) UTHIN(l) 


RUTH (NUTH) UTHIN(NUTH) 


Inlet Swirl Velocity 


NTT 


RTT(l) TTIN(l) 


RTT(NTT) TTIN(NTT) 

I NPT 

1 

j RPT(l) PTIN(l) 


RPT(NPT) PTIN(NPT) 


I UPMACH 


C Inlet Stagnation Temperature 


Inlet Stagnation Pressure 

or J characteristic values 
r 


166 


71 


76 


L 
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DT 


1 — r 


I 

*** UP TO 10 CARDS DESCRIBING 
*** OPERATOR SEQUENCE ARE NEXT. 

*** SEQUENCE IS TERLIINATED BY A BLANK 
*** CARD. SEQUENCE IS STORED IN NOR(I,J) 


lOP 


IML 


lEX 


UP TO 9 MORE DESCRIPTOR CARDS 


BLANK CARD 


NPl 


CODE 

JST 

JEND 

LST 

LEND 


***NP1-1 MORE CARDS DESCRIBING TYPE 
*** OF PRINTED OUTPUT DESIRED BEFORE THE 

NP2 

NTIM] 

SS 


CODE 

JST 

JEND 

LST 

LEND 



*** NP2-1 MORE CARDS DESCRIBING TYPE 

*** OF PRINTED OUTPUT DESIRED AFTER THE OPERATOR SEQUENCE 


NP3 


CODE 


JST 


JENDj 


LST 


LEND 


*** NP3-1 MORE CARDS DESCRIBING TYPE 

*** OF PRINTED OUTPUT DESIRED DURING THE OPERATOR SEQUENCE 
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Data Files Used for MESH3D, BLADE3D and GRAPH3D 

Several data files are used by MESH3D, BLADE3D and GRAPH3D as 
input files, scratch files, and permanent storage files. These files and 
their characteristics are shown below. 


Table 4.3. Data Files for MESH3D 


f"ortran 

Unit 

Number 

Type of File 

Comments 

' 

1 

Scratch 

UNFORMATTED 

Used by BLDT 

2 

Permanent 

UNFORMATTED 

Used by RGRID, INSEC, GETTHA, CALDR, CALDER 
for scratch file and BLDT for creating file 
TGEOM to be used by BLADE3D 

3 

Scratch 

UNFORMATTED 

Used by XRGRID, RGRID, BUILD 

4 

Permanent 

UNFORMATTED 

Used by GETTHA, CALDR, CALDER, 

BLDT as scratch file and BUILD for creating 
file GECM (to be used by MESH3D) 

5 

Input Data 
FORMATTED 

All Input Data in this file. Used by 
MESH3D, BLADE3D, XRGRID, FILL, RMID 

6 

Output Data 
FORMATTED 

Message and Error file used by all routines 

8 

Scratch 

UNFORMATTED 

Used by INSEC, GETTHA, CALDR, CALDER 

9 

Scratch 

UNFORMATTED 

Used by CALDER, BLDT 

12 

Permanent 

UNFORMATTED 

Used by SVGRID to create file SVSAVE which 
contains r,0,z coordinates of grid nodes 

IWR 

Output Data 
FORMATTED 

Used by all routines as an optional bulk 
output file 
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Table 4,4 Data Files for BLADE3D 


FORTRAN 

Unit 

Number 

Type of File 

Comments 

1 

Permanent 

UNFORMATTED 

Initial solution matrix storage 
used by START 

2 

Permanent 

UNFORMATTED 

Blade geometry for MESH3D, file TGEOM 
used by OPEN 

3 

Permanent 

UNFORMATTED 

Axisymmetric geometry file from MESH3D, 
file GEOM used by OPEN 

5 

Input Data 
FORMATTED 

All input data in this file, used by MAIN 
program and MTHREE 

6 

Output Data 
FORMATTED 

Used by all routines for Printed Output 

7 

Permanent 

UNFORMATTED 

Final solution matrix written by CLOSE 
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Table 4.5 Data Files for GRAPH3D 


FORTRAN 

Unit 

Number 

Type of File 

Comments 

1 

Permanent 

UNFORMATTED 

Solution Matrix File from BLADE3D 
unit number 7 

2 

Permanent 

UNFORMATTED 

File TGEOM from MESH3D unit number 2 

3 

Permanent 

UNFORMATTED 

File GEOM from MESH3D unit number 4 

4 

Permanent 

UNFORMATTED 

File SVSAVE from MESH3D unit number 12 

5 

Permanent 

FORMATTED 

Grid Numbering Information from BLADE3D 
unit number 5 

6 

Output Data 
FORMATTED 

Used for Printed Output 

7 

Permanent 

UNFORMATTED 

Solution Matrix Storage for Entire Finite 
Difference Grid 

8 

Permanent 

UNFORMATTED 

Solution Matrix Storage for S1-S2 Stream 
Surfaces 


Description of Printed Output for MESH3D 
The printed output from MESH3D consists of a message and error file 
which reproduces the axial grid numbering information, a running execution 
sequence log for major subroutines, and a bulk output file* This message 
file is of critical importance in locating input or logical errors in MESH3D 
operation. A sample message file is shown in Figure 24. The initial portion 
reproduces the axial grid input data, and if any parameter is outside of 
allowed ranges an error message will be printed in this file. Program execu- 
tion terminates when an input error is found. After successful completion of 
each principal subroutine a message is printed. This run file also contains 
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a set of warning messages from subroutine CALT which are generated when the 
blade row geometric input data must be extrapolated rather than interpolated. 
Printed values of blade row normal vectors and curvatures should be closely 
monitored to insure that the linear extrapolation performed represents the 
blade row geometry adequately. No further warning or error messages are nor- 
mally produced by CALT, BLDT or BUILD. 

The bulk output from MESH3D reproduces the geometric input data and 
the calculated results from MESH3D. Figure 25 illustrates the raw data input 
check for the hub and tip casing geometry, leading and trailing edge 
geometry, and damper geometry. This data is the "as read" geometry for each 
input data class. Figure 26 illustrates the scaled data smoothness check for 
XRGRID input. This data includes the scaled axial and radial coordinates and 
the first and second derivatives of radial position with respect to axial 
position. Small errors in these input quantities are usually quite apparent 
because of their influence on the second derivatives. The principal output 
from subroutine XRGRID is Illustrated in Figure 27 and consists of the axial 
and radial positions of the x and r coordinate line Intersections in physical 
space for each axial grid station. The calculated metric derivatives 
(3X/3z), (3X/3r), (3R/3r) and (3R/3z) are also shown. 

Bulk data output for subroutine CALT is the scaled blade row 
geometry as well as first and second derivatives of theta position with axial 
position for both pressure surface and suction surface. Best results from 
the calculation program BLADE3D are obtained when the first and second deri- 
vatives are smooth. When "jumps" in either derivative occur special care 
should be taken to Insure that the input data accurately represents the blade 
design intent. This bulk output is illustrated in Figure 28. The blade row 
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normal vectors and radii of curvatures calculated from the input geometric 
data are illustrated in Figure 29, This data is presented for each axial 
grid station inside the blade row. Special attention should be paid to the 
calculated radii of curvatures since small input errors greatly affect their 
value. This output is generated by subroutine BLDT. 

The final bulk, output from MESH3D is generated by subroutine BUILD 
and consists of the hub and tip slope and radius of curvature and the value 
of the transformation Jacobian, equation 23, for each axial and radial 
station. This output is illustrated in Figure 30. 

Printed Output for BLADE3D 

The printed output from BLADE3D begins by reproducing the input 
card file directing a particular solution pass. Figure 31 illustrates output 
from this first section. Particular attention should be paid to the printed 
MacCormack Operator sequence to insure that it is the one intended and that 
it is a symmetric sequence. At the user selected axial and radial node 
points, the solution matrix will be displayed for all theta nodes- Figure 32 
illustrates this output section. The content of this printed output, 
variables and coordinate system, is determined by options selected in the 
card input file for BLADE3D. If print option 4 or 5 is selected, then con- 
ventional theta averaged variables (meridional Mach numbers, total and 
pressure and temperature ratios, etc.) are calculated and displayed. This 
output is illustrated in Figure 33. Definitions of printed variables are 
given in special instruction number 4 for BLADE3D. 
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In addition to the user selected printed output, a substantial 
amount of information is printed during the operator sequence in order 
to document the history of a particular calculation. This information 
is printed after each 100 operator sequences and consists of printer/plotter 
plots of mass flow rate vs axial grid plane number, mid-span trailing edge 
static pressure vs operator sequence number and mid— span blade static 
pressure vs axial grid plane number. This information is the best guide 
to solution convergence rate and should be monitored closely. Examples 
of these three informational plots are shown in Figures 34, 35 and 36. 

Output Description for GRAFH3D 

The printed output from GRAPH3D consists of printer/plotter plots 
of blade surface MACH number and static pressure for each SI surface and 
plots of S2 streamsurf ace positions for the pressure surface, mid-channel 
surface and the suction surface. Examples of the SI and S2 surface plots 
are given in Figures 37 and 38, In addition, a summary of SI streamsurf ace- 
theta averaged flow variables is presented for the grid planes corresponding 
to the blade leading edge and the blade trailing edge. An example of this 
output is shown in Figure 39 , 

Both user defined solution matrix storage files have the same 
format which is specified in Table 4.6. Each FORTRAN record corresponds to 
one grid node point and contains the radial, tangential and axial grid index 
number, the non-dimensional radial, tangential and axial positions, the 
nondimens ional conservation variables (rp, rpu^, i^PUg, ^^t^ 

nondimensional physical flow variables (P , M , , M , T ), The finite 

difference grid output file contains all axial plane numbers 1 through NX, 

The S1-S2 streamsurf ace file contains only axial plane numbers N1 through N2. 
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EXAMPLE CALCULATION 

The computer codes described in the previous section were used to pre- 
dict the flow through a low aspect ratio, transonic compressor rotor de- 
signed by Urasek [22] of NASA LEWIS RESEARCH CENTER. This rotor has an 
inlet hub to tip ratio of 0,375, an aspect ratio of 1.56 and an inlet 
relative Mach number of 1.38, In conventional steady state testing at 
NASA LeRC, the rotor was found to have a peak adiabatic efficiency of 
90.6%, a total pressure ratio of 1.686 and a mass flow rate of 34.03 
Kg/sec. In Blowdown Tunnel testing at MIT [23], the rotor was found to 
have an adiabatic efficiency of 89.5%, a total pressure ratio of 1.677 
and a mass flow rate of 33,3 Kg/sec, 

The calculated operating point for the sample, inviscid calculation 
was set such that the predicted passage shock pattern approximated that de- 
termined in MIT BLOWDOWN COMPRESSOR TUNNEL flow visualization testing [24]. It 
was found to be impossible to match the mid-span bow-shock shape, the 
shape of tip wall static signature and the level of the tip wall static 
measurements. Predicted overall performance parameters were found to be 
a mass flow rate of 35.6 Kg/sec, a total pressure ratio of 1,756 and an 
adiabatic efficiency of 94,2%, Having generally the same shock structure 
at a higher mass flow rate and total pressure ratio is consistent with the 
assumption of inviscid flow. 

In order to verify that the calculated flow condition realistically 
represents the test rotor performance, comparisons between calculated and 
measured theta averaged performance is useful. A comparison of relative 
flow angle at the rotor trailing edge is shown in Figure 40 along with the 
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blade metal angles. For 0.90, the inviscid predictions follow the ex- 

pected variation with respect to the metal angles and match the experimental 
values. For R/R^^^>0.90, the experimental values show considerably less 
turning than would be expected either from deviation angle correlations 
or the inviscid predictions. The observed flow angle variation is 
difficult to explain on the basis of purely inviscid fluid phenomena, 
since as is shown in Figures 41 and 42, the experimental rotor leading 
edge meridional Mach numbers are nearly identical with the predicted 
values while the experimental total temperature ratios are larger than 
the predicted values. These phenomena suggest that either viscous 
phenomena such as profile boundary layer separation or viscous linked 
phenomena such as boundary layer flow migration significantly affect the 

flow for R/R ^>0,90. 
ref 

A comparison of the theta averaged rotor total pressure ratios is 
shown in Figure 43, The inviscid prediction and the NASA LeRC steady 
state measurements compare well over the entire span, while neither 
compares well with the high response total pressure measurements made 
in the MIT BLOWDOWN COMPRESSOR Tunnel . The origin of the high total 
pressure area, R/R^^^>0.70, in the MIT data, has not yet been explained. 

Measurements of the blade-to-blade variation in static density have 
been reported in reference [24], These studies have shown a rotor shock 
system consisting of a moderate strength bow shock at mid-span, a weak 
strength bow shock near the three-quarter span radius and of a secondary 
weak passage shock near the tip section. These characteristics are well 
illustrated in three contour plots, figures 44, 45 and 46, which show the 
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computed relative Mach number on blade-to-blade running surfaces. These 

surfaces have the constant radii of R/R „ = 0,80, 0.90 and 0.95. At an 

ref 

R/R^^^ of 0.80 a strong oblique bow shock is predicted with the flow remaining 
subsonic nearly everywhere downstream of the shock. At an R/R^^^ of 0,90 
a weak oblique bow shock is predicted, but the flow becomes supersonic in 
a small region near the trailing edge. At an R/R^^^ of 0,95, a weak bow 
shock is predicted with the flow remaining at supersonic across the entire 
passage. The supersonic region is terminated by a weak compression or shock 
and the flow is subsonic everywhere downstream of the trailing edge. 

To clarify the predicted shock structure, relative Mach number contour 
plots for the pressure and suction surface are shown in figures 47 and 48, 
Figure 47 shows the flow to be subsonic over nearly the entire span. The 
supersonic flow portion is terminated by a weak shock. The suction surface 
flow pattern is much more complex as shown in figure 48, The flow is super- 
sonic over a large fraction of the blade surface with two different passage 
shocks evident. The suction surface shock structure is sketched in figure 
49 which illustrates that the first passage shock, really the bow shock, 
is a strong oblique shock near mid-span and a weak oblique shock near the 
tip section. The two shock families merge near the three-quarter span point 
where the shock turning angles are the same for both families. 
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The intra-blade flow visualization technique does not allow resolution 
of leading edge flow features for R/R^^^>0,85 and high response wall 

static measurements must be used to experimentally probe the shock structure 
near the tip radius. Reference [23] presents wall static measurements at 
the axial locations illustrated in Figure 50, and sample static pressure 
traces are reproduced in Figures 51, 52 and 53, The predicted wall static 
pressures are also compared to 5 cycle ensemble averages of these data in 
Figures 54, 55 and 56. At the upstream measurement ports 2.0 and 56, the 
shape of the computed wall static pressure trace closely follows the measured 
5 cycle average, but the shock pressure rise is well under-predicted. At 
the mid-chord measurement port 2,5, the shape of the pressure trace is well 
predicted, but the mid-passage pressure level is predicted to be too high. 

In order to clarify the comparison of predicted and measured wall 
static pressures, the port 56 measurement is compared to the computed leading 
edge pressure in Figure 57, This figure shows that the correct bow shock 
pressure rise is predicted, but the predicted bow shock appears at the wrong 
axial position. Since the predicted shock pressure rise is consistent with 
the blade leading edge wedge angle, it must be concluded that the tip bow 
shock is detached in the experiment. The most likely explanation for this 
difference is that the experimental compressor has a tip "end-bend" or 
local over-twist to accommodate a tip end-wall boundary layer. The end- 
wall boundary layer is absent in the MIT test configuration. 
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To complete the documentation of the sample solution, blade surface 
static pressures and blade surface Mach number distributions along each 
radial grid plane are shown in Figures 58 through 87, Figures 58 
through 72 illustrate the static pressure, and Figures 73 through 87 
illustrate the relative coordinate system total Mach number. The nominal 
upstream sonic line occurs at radial plane number 9, but supersonic flow 
extends inward to radial plane number 5. These figures illustrate the 
transition from subsonic flow to supersonic flow; the inviscid calculation 
leading edge resolution in subsonic flow. Figures 59 and 73; and the shock 
resolution. Figure 87. These figures show that the compromises in 
leading edge resolution have not greatly degraded the solution. 



Table 4.6, Format of Solution Matrix Storage File Produced by GRAPH3D 


RECORD 

NUMBER 

TANGENTIAL STATION 
NUMBERS 

3 INTEGER VALUES 
£ k i 

RADIAL, TANGENTIAL, 
AXIAL POSITIONS 
3 FLOATING POINT 
VALUES 

FLOW VARIABLES 

rp, rpu , rpu„, rpu , rpE 
r U z 

5 FLOATING POINT VALUES 

FLOW VARIABLES, LABORATORY 
COORDINATES Mg, , T^ 

5 FLOATING POINT VALUES 

1 

1 

1 

1 

0.3093, 0.,-C.4880 

0.270, 0., 0. , 0.138, 0.493 

0.714, 

0., 0., 0.547, 1.00 

2 

1 

2 

1 

0.3093, 0., -0.4880 

same 


same 

NTH 

1 

NTH 

1 

0.3093, 0., -0.4880 

1 

same 

1 

same 

NTH + 1 

2 

1 

1 

0.3904, 0., -0,4790 

0.337, 0. , 0. , 0.179, 0.614 

0.714, 

0., 0., 0,547, 1.0 

2 NTH 

2 

NTH 

1 

0.3904,0., -0.4790 

same 


same 

NTH * R 

NR 

NTH 

1 

1.0407, 0., -0.4180 

0.909, 0., 0. , 0.466, 1.658 

0.714, 

0., 0., 0.539, 0.99 

1 + NTH * R 

1 

1 

2 

0.3801, 0., -0.4337 

0.379,-0.10, 0.0,0.18, 0.519 

0.714, 

-0.035, 0., 0.612, 1.0 

2 * NTH * NR 

NR 

NTH 

2 

1,0408, 0. , -0.3536 

0.817,-0.006, 0.001,0.442, 1.55 

0.714, 

-0.015, 0.0, 0.546, 1.0 

NX * NTH * NR 

NR 

NTH 

NX 

.9796, 0.0, 0.8463 

1,231, 0,, 0.426, 0,631, 2.645 

1.101, 

0,0, 0.345, 0.479, 1.150 
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FIGURE 1 


PHYSICAL SPACE VIEW OF SINGLE BLADE PASSAGE 
SHOWING COMPUTATIONAL DOMAIN 
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R = f ^ 0.5 

“"tip ^hub 2 te ■" ^ A 


FIGURE a 


PLAN VIEW OF PHYSICAL SPACE AND 
FIRST COMPUTATIONAL DOMAIN 
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FIGURE 3 


BLADE TO BLADE VIEW OF PHYSICAL SPACE 
AND FIRST COMPUTATIONAL DOMAIN 
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FIGURE 4 


PLfiN VIEW OF ROTOR GEOMETRY SHOWING 
PflRT-SPfiN SHROUD 
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FIGURE 5 PLRN VIEW OF COMPUTRT I ONflL DGMfilN 
SHOWING PART-SPAN SHROUD 



RADIUS 

~ 0 .< t 0 - 0.20 0.00 0.20 0 .^Q 0.60 0.80 1.00 


GRID STRUCTURE IN THE R - Z PLfiNE 



FIGURE 6 
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FIGURE 7 


GRID STRUCTURE IN R - THETR PLANE 
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k = NTH 
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Blade Pressure 
Surface 
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Periodicity Enforced 


FIGURE 8 


BLRDE TO BLADE VIEW OF COHPUTflT lONflL DOMfilN 
SHOWING DIFFERENT BOUNDRRT TTPES 
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FIGUBE 9 


EXAMPLE GRID FOR CONTINUITY INTEGRATION 
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SURFfiCE NORMAL VECTOR 

UNIT VECTOR TANGENT TO 3 COORDINATE LINE 
UNIT VECTOR TANGENT TO />j COORDINATE LINE 


FIGURE 10 COORDINATE SYSTEM FOR WALL STATIC 
PRESSURE EVALUATION 




|(i + l,3) 


Solid 

Wall 



FIGURE 12 GRID POINT NUMBERING SCHEME FOR STAGGERED GRID SYSTEM 



Dummy Grid Line 
K = 1 


Periodic — v 
Boundary \ 

W^J(\ 


Blade Surface 


Periodic 

Boundary 


\ /YY Y\\VCYy \Y 

Ay^\ \ \'\ y0"^^\ 


/I U l> i 




Dummy Grid Line 
K = NTH 


Blade Surface 


/f / i/'^r 

\yXy 


\ — Periodic Boundary 

FIGURE 13 BLfiOE TO SLADE VIEW OF STAGGERED GRID SYSTEM 




2 3 Direction 


FIGURE m INFLOW BOUNDARY CHARACTERISTIC FORMULATION 
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END 

FIGURE 17 FLOWCHART FOR PROGRAM MESH30 
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FIGURE 18 FLOWCHART FOR PROGRAM 8LRDE3D 









Trai ling 
Edge 


Finite Difference Grid Lines 

Blade Specification Lines 

Note ' Each Blade Line Must Have The 
Some Number of Input Values. 


FIGURE 19 PLAN VIEH OF BLfiDE GEOMETRY SHOWING TYPICAL POSITIONS 

OF FINITE DIFFERENCE GRID LINES AND BLADE SPECIFICATION 



RADIUS 
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KX=0.001 



-0.50 -0.30 -0.10 0.10 0.30 

AXIAL DISTANCE 
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0.70 


0.90 


FIGURE 20 flXIRL GRID PLANE SPACING EXAMPLE 
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FIGURE 23 


BLRDE GEOMETRY INPUT TRBLE fJEF I N I T I ONS 



7jI\ mXJAL OtiMPRtSSOR OHlh CDDR ♦ 

-NASA I..OW ASPECT rat;i:o GK;:n;t test 

-THE NUMBER OF BLADES IS 22 

-HALF MESH SPACING (DjNSTANf IS 0.000 

-GRin POINTS IN X ):HKE(;:T;i ON 30 

-IN THETA DIRECTION ■■: 10 IN R iHPEiCnON 10 

-THE BI.AHF LEADINC EDCI: AND TRAILING EDGE ARE AT AXIAI.. CiRID LINES 3 0 20 

-KUTTA point is at liRHi t.INE 20 

-THE FIRST POINT ON THE DAMPER IS 0 

-THE LAST POINT ON THE DAMPER IS 0 

"OU’IPin KILE NUMBER IS 'lA 

-DEBUG DUMP PARAMETER IS 0 


SUBROUTINE GRID FINISHED. 
EXTRAPOLATION NEEDED FOR 

IER^^-= 

1. 

0 
;i 0 

1 

EX TR APT.)!.. AT ION 

NEEDED FOR 

■J 7 1.. 

u 

1 

EXTRAPOLATION 

NEEDED FOR 

,.l 7 1. 

12 

1 

EXTRAPOLATION 

NEEDED Ef)R 

■J 7 L 

1.3 

1 

EXTRAPOLATJ ON 

NEEDED FOR 

J 7 L. 

14 

1 

EXTRAPOLATION 

N E E it E Li E 0 R 

J 7 )._ 

ILi 

1 

•SUBROUTINE CAL 

.1 FINISHED. 

lER^^ 

0 



-SUBROUTINE BLDT FINISHED 
-SUBROUTINE BUILD FINISHED 
"SUBROUTINE SOGRID FINISHED 
-STOP 


FIGURE 34 MESH3D RUN LOG FILE 
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HUB AND TIP COORDINATE INPUT CHECK 


HUB? TIP? l > e *7 A?n;i r»t:> cooKiHiiATES » 


SCftLt.D RAW DAfA INPUT XrY PAIRS 


>12771 

0 > 3!:iOG2 

-1.3 7908 

0.3[=i0l;i2 

-0 >ri:>oi / 

0 . 3^^052 


0.34697 

-0 , \ 2A:^2 

0 » 3 U 3 3 

-o.or/Bi 

0,37192 

0 • f ) / \'i ^ O 

0 > 3 9 U 9 6 

0 . J 71-?65 

0,42765 

0 . 27 J.64 

0.45072 

0 > 7 0 6 3 

0*46335 

0. ^6962 

0.47309 

0 . L^i6B61 

0.48474 

o.r;/4Av 

0. 13474 

0.613AS 

0.43474 

0 . 6:i26 4 

0.43174 

0 . 3 <^>1 

0.43474 

1 . I^i928 

0.43174 


3CAI..H:Li raw data input x?y pairs 


-3. 3.2332 

1. , 00000 

- 1 . J. 7 9 3 1 

1 .00000 

"0.52027 

1.00000 

"0 . 32237 

1.00000 

-0 .12435 

1 . 00000 

-0.01781 

0.999B8 

0 . 07367 

0.99673 

0 . :l 7263 

0.97950 

0.27169 

0.96207 

0. 37070 

0.95346 

0 . 46971. 

0.95049 

0 . 56872 

0.95049 

0 ♦ 57430 

0.95049 

0 .61 378 

0.95049 

0 . 65276 

0.95049 

0 . 691. 7 4 

0.95049 

1 .. J. 595 1, 

0. 95049 


FIGURE 25 


MESH3D BULK OUTPUT FILE 


HUB AND TIP COORDINATE SMOOTHNESS CHECK 


HU£tp 

ilh'f L.*lr:.r 

AND T.t:. 

I'lERIOIT lots. 


J 

X 

Y 

OYDX 

D2YDX 

1 

~3« 1283 

1.0000 

0.0000 

0.0000 

'y 

im 

“1 » 1793 

1.0000 

0.0000 

0.0001 

:i 

-0.5203 

1.0000 

-0.0001 

-0.0005 

A 

-0,3224 

1.0000 

0.0003 

0.0041 

5 

-0.1243 

1.0000 

-0.0009 

-0.0159 

6 

-0,0178 

0,9999 

-0.0006 

0.0215 

7 

0.0737 

0.9967 

-0.1033 

-2.2658 

8 

0 , 1727 

0,9795 

-0.2033 

0.2448 

9 

0.2717 

0.9621 

-0.3334 

3.3673 

10 

0 , 3707 

0.9535 

-0.0520 

0.4784 

11 

0.4697 

0.9505 

- 0.0095 

0.3784 

12 

0,5687 

0.9505 

0.0003 

-0.1786 

13 

0.574S 

0.9505 

-0.0001 

0.0136 

14 

0.6138 

0 .9505 

0.0000 

“0.0036 

15 

0.6528 

0.9505 

0.0000 

0.0009 

16 

0.6917 

0.9505 

0.0000 

0.0000 

17 

1 .1595 

0.9505 

0.0000 

0.0000 


FIGURE 26 MESH3D BULK OUTPUT FILE 



R-2 COORDINATES IN PHYSICAL SPACE 


X-R CQORiUMHTtiH AN);i Mi:: i kites, 

THE AXIAL PACKIWe FAClOR IS 0,100 


AXIAL 

RAKaAI.. 

AXIAL 

RAHIAL 

STATION 

S i A 1 1 ON 

POSITION 

POSil CION 

1 

1 

-0, 49A4 

0,3094 

;L 


-0,4876 

0, 390U 

1 

s 

-o,48;r/ 

0 , A 7 1 7 

1 

A 

-0,4776 

0,:;)029 

;i 

I'i 

-0,4690 

} 6 3 4 J. 

1 

6 

-0,4U99 

0,7133 

1 

7 

-0, 4r^2:i. 

0 , 7 9 6 6 

1 

8 

-0,4460 

0,8779 

1 

9 

-0, !140:I. 

0 , 9 5 9 3 

1 

10 

-0,4::; 78 

.1 ,0407 


XX 

Xk 

Kk 

RX 

3,538 

— 0 ,436 

1,339 

0 ♦ 0 3 3 

3,4 97 

-0,224 

1 , 338 

0.0::i9 

3,463 

-0 ,209 

1 ,533 

0,026 

3 ,416 

-0*289 

1 , 5 3 S 

0,022 

3, 3!:; 7 

“0,428 

J. ,537 

0,018 

3,292 

-0,327 

1*537 

0,014 

3 , 2 3 4 

-0,267 

1 , 2* 3 6 

0,010 

3< :l 89 

-0,203 

J * 1 1 3 6 

0,006 

3 , 1 8 

-0*42 0 

1 * 5 3 6 

0,002 

3,041 

— i'\ A ’J L 
* V u 

■1 r* -T r* 
1 ^ 

-0,002 


FIGURE 27 


HESH3D BULK OUTPUT FILE 
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no 

SCALED BLADE ROW CEOMERTY 

NASA LOW ASPLOT RATIO TWO-STAOE -- NEW ROTOR ONE GEOMETRY — 8/23/73 
INPUT DATA FOR SEOTION 1 


L 

J 

X 

R 

RTHPS 

RTHSS 

1 

1 

0.00000 

0.40147 

0.00160 

-0.00160 

1 

2 

0. 00127 

0. 40126 

0.00075 

-0.00315 

i 

3 

0.03 699 

0.39901 

-0.00994 

-0.02000 

1 

4 

0,03271 

0.39743 

-0.01973 

-0.03438 

1 

S 

0.04843 

0.39634 

-0.02836 

-0.04727 

1 

6 

0*0A4i:.> 

0.39563 

-0.03603 

-0.05870 

.1 

7 

0.07938 

0.39519 

-0.04301 

-0.06863 

1 

a 

0. 09T>A0 

0 , 39496 

-0.04931 

-0.07719 

1 

9 

0.11332 

0.39437 

-0.05500 

-0.08454 

1 

10 

0.12704 

0,39433 

-0.06003 

-0.09073 

1 

] 1 

0.14276 

0.39496 

-0.06442 

-0.09583 

1 

12 

0,irJ843 

0.39505 

-0.06813 

-0.09984 

1 

13 

0.3 7420 

0.39516 

-0.07115 

-0.10279 

1 

14 

0 . 13993 

0.39524 

-0.07343 

-0.10468 

1 

10 

0.20565 

0.39529 

-0.07495 

-0.10546 

1 

1 ^ 

0.22137 

0.39530 

-0.07566 

-0.10512 

1 

1 7 

0.23709 

0.39526 

-0.07550 

-0.1 0355 

1 

13 

0.25231 

0.39513 

-0.07442 

-0.10069 

1 

19 

0 . 26853 

0.39506 

-0.07233 

-0.09641 

1 

20 

0.28426 

0.39494 

-0.06912 

-0.09057 

1 

21 

0.29998 

0.39487 

-0.06469 

-0.08293 

1 

2 2 

0.31570 

0.39491 

-0.05334 

-0.07315 

1 


0.33142 

0.39518 

-0.05135 

-0.06085 

1 

24 

0.34714 

0.39539 

-0.04230 

-0.04448 

.1 

25 

0.34805 

0.39595 

-0.04179 

-0.04331 


SECOND DERI DATIVE CHECK FOR PRESSURE SURFACE 


J 

X 

Y 

DYDX 

D2YDX 

1 

0.0000 

0.0016 

-0.6722 

-1.7001 

2 

0.0013 

0.0008 

-0.6754 

-3.4001 

3 

0,0170 

-0.0099 

-0.6619 

5.1238 

4 

0.0327 

-0.0197 

-0.5845 

4.7202 

5 

0 >0434 

-0.0284 

—0.5155 

4.0563 

6 

0.0642 

-0.0360 

-0.4644 

2.4397 

7 

0.0799 

-0.0430 

-0.4224 

2.9042 

8 

0.0956 

-0.0493 

-0.3810 

2.3609 

9 

0 ,111 3 

-0.0550 

-0.3412 

2. 7080 

10 

0.1270 

-0.0600 

-0.2998 

2.5585 

11 

0.1428 

-0.0644 

-0.2579 

2.7705 

12 

0.1585 

-0.0681 

-0.2141 

2.7992 

13 

0 . 1742 

-0.0711 

-0. 1637 

2.9773 

14 

0.1899 

-0.0734 

-0.1212 

3.0727 

15 

0 . 2056 

-0.0749 

-0 . 0716 

3.2260 

16 

0.2214 

-0.0757 

-0.0182 

3.5723 

17 

0,2371 

— 0 . 0755 

0.0388 

3.6729 

18 

0.2528 

-0.0744 

0.0996 

4.0731 

19 

0.2635 

-0.0723 

0.1674 

4.5474 

20 

0 . 2843 

-0.0691 

0.2414 

4.8666 

21 

0.3000 

-0,0647 

0,3252 

5.7928 

22 

0.3157 

-0.0588 

0.4195 

6.2089 


FIGURE 28 


MESH3D BULK OUTPUT FILE 



CflLCULRTED BLADE NORMAL VECTORS AND CURVATURES 


NORMAL UEC' 

roRS AND rad:i:i 

UF CURVATURE 

FUR a>::lal 

STATION 

K'O, 

10 


L COSNR 

COSNT 


RES 

RTAU 

1., COSNR 

COSNT 

COSH^. 

R E S 

RTAU 

PRESSURE SIDE NORMAL VECTOR 







1 “0,099 

0.841 

0,532 

0,000 

0,000 

2 -0,073 

0,823 

0< 563 

0,000 

-0,198 

4 -0,1.48 

0,765 

0,627 

0,000 

0,074 

5 -0,142 

0,678 

0,721 

0 , 000 

0,350 

7 -0.070 

0,573 

0,817 

o 

o 

o 

0,299 

8 -0,060 

() ,532 

0.845 

0,000 

0,686 

10 -0,072 

0,435 

0,897 

0,000 

0,000 








FIGURE 

29 

MESH3D 

BULK OUTPUT FILE 





CALCULATED HUB AND TIP CASING SPLOE 
CALCULATED TRANSFORM JACOBIAN VALUES 


HUB 8L0PL HUB CURV* 
-0,02018 0,04381 

JACOB IANS FOR RADIAL. GRID 
;l ;l9»22Ji62 :> :18,8J300 

7 17,41035 8 17,15939 


TIP S1..DPF 
-0,00006 
OCAT J ONS, 

3 18,66972 

9 16,88091 


•HP CURV, FOR A>:iAL STATION NO, 
0,00151 

4 18,41212 5 18.09552 6 17,7 

10 16,34608 


FIGURE 30 


MESH3D BULK OUTPUT FILE 


1 

r!968 
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112 


INyiSCID 3-H AXIAL COMPRESSOR ANALYSIS CODE* 
NASA LOW ASPECT RATIO ROTOR TEST 
THE NUMBER OF BLADES IS 22 
50 cycles are to be run with DT --- 0.00200 


OPERATOR SEQUENCE IS 

OPERA f OR tDT MULT IPLERf repeat COUNT 


r 

DP 



1 


8 

X 

np 



1 


a 

R- 

0 F^' 



8 


1 

R 

0 i ■’ 



8 


1 

X 

OP 



1 


8 

T 

OP 



1 


8 

GRID 

POINTS 

IN 

X 

DIRECTION 

IS 

60 

GRID 

POINTS 

IN 

THETA DIRECTION 

IS 

G R“ .L I-i 

POINTS 

IN 

R 

DIRECTION 

IS 

18 


THE BLADE LEADING AND TRAILING EDGES ARE AT AXIAL STATIONS 25 AND 45 

the: KUTTA point is located at axial station 45 

FIRST DAMPER STATION IS AT 0 

THE LAST DAMPER STATION IS AT 0 

lULOW IS 1 

IDHIGH IS 18 

GAMMA IS 1.400 

BLADE SPEED (WRT/AO) IS 1.271 

ARTIFICAL UISCOSITY PARAMETER IS 0.400 

DOWNSTREAM PRESSURE IS TO BE SET TO 0.780 

UPSTREAM MACH NUMBER PARAMETER IS 0.510 


L 

RADIUS 

UTHETA 

1 

0.33194 

0,00000 

**l'I 

0*37240 

0,00000 


0.41238 

0.00000 

4 

0.45336 

0.00000 

5 

0.49385 

0.00000 

6 

0 . 53432 

0 . 00000 

7 

0.57480 

0.00000 

8 

0 . 61528 

0 ,00000 

9 

0.65576 

0 . 00000 

10 

0.69625 

0.00000 

1 1 

0.73674 

0.00000 

12 

0 . 77723 

0.00000 

13 

0.31773 

0.00000 

14 

0 . 85823 

0.00000 

1 5 

0.89874 

0.00000 

16 

0 i 9 3 9 2 4 

0 . 00000 

17 

0 . 97975 

0.00000 

18 

1 ♦ 0 2 0 2 6 

0.00000 


TTIN PTIN OR JPLUS 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5,32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5,32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5.32000 

1.00000 5,32000 


FIGURE 31 


BLRDE3D RUN LOG FILE 



GRID (3UTPUT FOR 30 TO 30 IN APS SYSTEM 
NASA LOW ASPECT RATIO ROTOR TEST 


30 5 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

PT 

0*944 

0*940 

0*937 

0*934 

0.932 

0*929 

0*925 

0*919 

0*912 

0.904 

MR 

0*102 

0*102 

0.103 

0.105 

0*106 

0.108 

0.112 

0.115 

0 * 1 1 7 

0.119 

MTH 

0*355 

0*357 

0*351 

0*345 

0.33 7 

0*327 

0.31 6 

0 . 3 0 5 

0.294 

0.282 

MX 

0*521 

0*523 

0*533 

0*547 

0.561 

0.576 

0 ♦ 539 

0.60 2 

0.614 

0*625 

TT 

1*116 

1*109 

1*103 

1.097 

1.091 

1 *085 

1*081 

1*078 

1*075 

1*072 

THETA AVERAGE PROPERTIES 

♦ 








XSPAN 

PT 

MR 

MTH 

MX 

TT 






0,227 

0*893 

0*119 

0*294 

0 i 6 0 

1.07 4 







FIGUBE 32 BLBDE3D BULK OUTPUT FILE 

SOLUTION MfiTRIX FOR 

fiXIflL GRID PLfiNE 30 
RADIAL GRID PLANE E 
NON-DIMENSIONAL FLOW VARIABLES AT 
TflNGENTRL POINTS 1 TO 10 



114 


[._ 

Radius 

H A C l-l 

U M B E. R S 


P ON 

T ON 

ANGLES 




MERJ.D 

TANQEN 

TOT AL 

PREF 

TREE 

WHIRL 

HER ID 

:i. 

0*381 

0*556 

0*34 

0*654 

1*257 

1 * 002 

1 ♦ / 

18*1 


0*419 

0*544 

0*31 

0*624 

1 , 199 

1*060 

29 * 4 

16*8 

,5 

O' * 4 3 

0*577 

0*31 

0*6 5 6 

1 + 222 

1 ♦ 062 

27*9 

14 * 3 

4 

0 * 4 V 6 

0*599 

0 * 30 

0*671 

1 * 2 4 1 

1*069 

26*9 

1 2 * 6 

1L“ 

0 * 3 3 3 

0*615 

0 * 29 

0 * 682 

1*259 

1 * 0 7 4 

25*5 

11*1 

W 

0*571 

0*625 

0 ♦ 2 Q 

0*685 

1*273 

1 * 0 7 7 

24*2 

10*1 


0 * 608 

0*630 

0*26 

0*683 

1*280 

1*077 

22*7 

9*4 

P, 

0*645 

0*635 

0v24 

0*673 

1*269 

1*073 

20*5 

3*5 

9 

0* * o c) .1. 

0 * 647 

0*20 

0 *673 

1*245 

1*065 

1 7 * 6 

7*0 

1 0 

0*718 

0*660 

0*17 

0*681 

1*215 

1*057 

14*2 

4 * 8 

i 1 

0*733 

0*66 3 

0*14 

0*679 

1 * 185 

1 * 049 

11*5 

2 * 8 

12 

0*791 

0 ♦ 6 6 

0*11 

0*673 

1*157 

1*042 

9 * 4 

1 * 0 

.1 O 

0*827 

0 * 65 9 

0*09 

0*666 

1 * 134 

1*035 

7*9 

"0*8 

14 

0*363 

0*652 

0*08 

0*657 

1*114 

1*030 

6*7 


13 

0*399 

0*640 

0*07 

0*643 

1*098 

1*027 

5 * 8 

“4*7 

16 

0‘* 935 

0*627 

0*06 

0*629 

1*078 

1*025 

5*1 

- 7 *5 

17 

0*970 

0*597 

0*06 

0*600 

1 * 101 

1*024 

5*4 

“7*4 

1 B 

1*005 

0*589 

0*06 

0*592 

1*066 

1*034 

5 * 7 

■■ 1 w * 6 


FIGURE 33 BLRDE3D BULK OUTPUT FILE 

THETfl flVERflGED FLOW QURNTITES 
RXIRL GRID PLRNE 30 



AXIAL GRIli PL. ANE 


:i *ooE+oo 


•j. >08 h:+oi 

2*07E+01 

;s*ouE-i-oi 

— — — — — X — — — 
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